Residual Decomposition for Lithotype-Aware Characterization of Rock Mechanical Parameters from Well Logs Under Lithological Heterogeneity

1. Introduction

As energy exploration advances into deeper coal-bearing formations, ensuring wellbore stability and operational safety has become increasingly challenging [1,2,3]. In such environments, strong lithological heterogeneity, complex cleat systems, and frequent coal–rock interbedding lead to substantial variability in geomechanical properties. Reliable characterization of parameters such as uniaxial compressive strength (UCS), elastic modulus, and Poisson’s ratio is therefore essential for drilling design, wellbore stability evaluation, and hydraulic fracturing optimization. However, in practice, these parameters remain highly uncertain, significantly increasing operational risks and compromising engineering reliability [4,5,6].

The core difficulty of this problem is not merely nonlinearity, but non-uniqueness induced by heterogeneity. In heterogeneous formations, similar well logging responses may correspond to fundamentally different mechanical properties due to variations in lithological composition and structure. This leads to an ill-posed inverse problem, where a single global mapping from logging data to geomechanical parameters does not exist. As a result, any method that assumes a unified mapping inevitably introduces systematic bias.

Existing approaches for rock mechanical parameter characterization can be broadly categorized into three groups. First, empirical methods establish predefined relationships between logging responses and mechanical properties [7,8]. While computationally efficient, they rely on strong assumptions of formation homogeneity and often fail in heterogeneous intervals. Second, data-driven methods, including machine learning models such as random forests, gradient boosting, and neural networks, aim to capture complex nonlinear relationships [9,10]. Third, feature-augmented learning approaches incorporate additional geological information, such as lithology indicators, as input features to improve prediction accuracy.

Despite their differences, all these approaches share a fundamental limitation: they implicitly assume a single global mapping between input features and mechanical properties. This assumption is incompatible with heterogeneous geological systems, where lithological regimes introduce conditional deviations in geomechanical response. Consequently, even advanced machine learning models tend to average over lithological differences, leading to biased predictions, particularly in lithological transition zones.

From a geological and engineering perspective, lithotype variations are a primary factor controlling mechanical behavior in coal-bearing formations [11,12,13,14]. Differences in maceral composition, fracture density, and structural integrity result in distinct strength and deformation characteristics across lithotypes. Ignoring these variations not only reduces predictive accuracy but also violates physical consistency, as the model fails to reflect the underlying geological controls on mechanical response.

To address this fundamental limitation, we propose a heterogeneity-aware modeling framework based on lithotype-conditioned residual learning. Instead of treating lithotype as an auxiliary feature, the proposed method formulates geomechanical characterization as a structured modeling approach, in which the overall response is separated into a global component and a lithotype-dependent residual. This formulation explicitly models conditional deviations induced by lithological regimes, thereby helping address heterogeneity-induced non-uniqueness. This is not a feature augmentation strategy, but a structured modeling of lithotype-dependent bias.

The proposed framework is potentially extendable to other heterogeneous subsurface systems and applicable to a wide range of heterogeneous geological systems where conditional bias exists. By explicitly modeling lithotype-dependent deviations, it restores physical consistency in prediction and mitigates systematic bias introduced by global mapping assumptions.

The effectiveness of the proposed method is evaluated through cross-well validation and comparative analysis. The results demonstrate that the method consistently reduces lithotype-induced bias and achieves stable generalization across different wells and geological conditions. Further analysis shows that the residual component exhibits structured patterns aligned with lithological regimes, confirming that the performance gain arises from modeling conditional bias rather than introducing additional information.

The main contributions of this study are summarized as follows: (1) A heterogeneity-aware residual learning paradigm is proposed to address non-uniqueness in geomechanical characterization under lithological heterogeneity; (2) A structured decomposition formulation is introduced to explicitly model lithotype-induced conditional bias; (3) The proposed method demonstrates improved cross-well generalization by resolving heterogeneity-induced prediction bias; (4) The physical consistency of the model is validated by linking residual behavior to lithological regimes.

Overall, this study reframes rock mechanical parameter characterization as a structured learning problem under heterogeneity, providing a generalizable and physically consistent modeling framework for complex subsurface systems.

2. Materials and Methods

This study proposes a heterogeneity-aware modeling framework for rock mechanical parameter characterization under strong lithological variability. Unlike conventional workflows that rely on feature integration, the proposed approach formulates the problem as a structured modeling approach, in which geomechanical response is decomposed into a global component and a lithotype-conditioned residual. This formulation enables explicit modeling of heterogeneity-induced conditional bias rather than implicitly absorbing it into a single global mapping.

Specifically, lithotype information derived from the HMLZ index is not treated as an auxiliary feature, but as a lithotype variable that governs systematic deviations in mechanical response. By integrating logging features with lithotype-conditioned residual modeling, the method captures cases where similar logging signatures correspond to different mechanical properties under distinct lithological regimes. As such, the proposed approach constitutes a model decomposition of geomechanical response under heterogeneity.

The overall methodology consists of three components. First, lithological heterogeneity is quantified through lithotype identification based on the HMLZ index, providing a structured representation of geological regimes along the wellbore. Second, a global mapping is learned from logging responses to geomechanical parameters, capturing lithology-independent trends. Third, a residual model conditioned on lithotype is constructed to learn systematic deviations induced by heterogeneous lithological conditions. The final prediction is obtained by combining the global and residual components, yielding a physically consistent characterization of rock mechanical behavior.

This formulation is potentially extendable to other heterogeneous subsurface systems and can be applied to other heterogeneous geological systems with similar characteristics, extending beyond coal-bearing formations.

2.1. Study Area Overview

The study area is located in the Ordos Basin, a large cratonic sedimentary basin in the western part of the North China Platform, with a sedimentary thickness of approximately 5,000 m. The Upper Paleozoic strata are dominated by fluvial–deltaic clastic deposits and host extensive coal-bearing formations [15]. The investigated region lies within the Daniudi Gas Field, a key area for deep coalbed methane (CBM) exploration and development.

Coal seams in this region are typically buried at depths exceeding 2,000 m and exhibit pronounced lithological heterogeneity characterized by frequent coal–rock interbedding and well-developed cleat systems. As drilling activities extend into deeper formations, engineering challenges such as wellbore instability and borehole enlargement have become increasingly severe. These challenges are closely associated with lithotype-dependent variations in geomechanical properties, which introduce significant uncertainty into stability evaluation and drilling design.

Within this geological context, the characterization problem is inherently heterogeneous, where identical or similar logging responses may correspond to different mechanical properties across lithological regimes. This makes the study area a representative testbed for evaluating methods designed to help address heterogeneity-induced non-uniqueness.

Accordingly, this study adopts the Ordos Basin dataset to validate the proposed heterogeneity-aware modeling framework, with a focus on assessing its ability to capture lithotype-controlled deviations and to achieve reliable geomechanical characterization under complex geological conditions.

2.2. Dataset Construction

To support the proposed heterogeneity-aware modeling framework, a cross-well dataset was constructed to explicitly preserve lithotype-induced variability and enable the modeling of lithotype-dependent residuals under heterogeneous geological conditions. The dataset is designed not merely for predictive accuracy, but for capturing the non-uniqueness of the mapping between logging responses and geomechanical properties.

Each depth sampling point is treated as a fundamental unit of analysis. For a given depth d, the corresponding logging responses are represented as the input vector $X\left(d\right)$, while the measured geomechanical parameters form the target vector $Y\left(d\right)$. In addition, lithotype information $L\left(d\right)$ derived from the HMLZ index is incorporated as a lithotype variable. The resulting data structure $\left(X\right(d),L(d),Y(d\left)\right)$ enables explicit modeling of lithotype-dependent deviations, which is essential for resolving heterogeneity-induced ambiguity.

The input features consist of conventional well logging curves, including acoustic transit time (AC), bulk density (DEN), gamma ray (GR), deep and shallow lateral resistivity (LLD, LLS), and spontaneous potential (SP). These measurements collectively characterize formation properties from complementary physical perspectives. Importantly, lithotype is not treated as a simple auxiliary feature, but as a variable that governs conditional shifts in the mapping from X to Y. This distinction is critical: the objective is not feature enrichment, but the representation of structured bias induced by heterogeneous lithological regimes.

The definitions of input features and output targets are summarized in Table 1.

The target vector is defined as a five-dimensional geomechanical parameter set:

$$Y\left(d\right)=\left[E_s\left(d\right),\nu \left(d\right),UCS\left(d\right),C\left(d\right),\varphi \left(d\right)\right]$$

(1)

This multi-output formulation preserves intrinsic correlations among mechanical parameters and allows consistent characterization across depth. Crucially, under heterogeneous conditions, the mapping from X to Y is not unique; instead, lithotype L induces structured deviations that must be explicitly modeled.

Data preprocessing ensures consistency and reliability of the dataset. All logging curves are aligned in the depth domain through datum unification and sampling interval standardization. Quality control procedures include the identification of missing intervals, suppression of spike noise, and removal of physically implausible values. Minor gaps are filled via interpolation to maintain continuity without altering global trends. After preprocessing, all features and targets are mapped onto a unified depth grid, forming a cross-well consistent data matrix.

A well-based partitioning strategy is adopted to ensure rigorous evaluation of generalization capability. By treating each well as an independent unit, the dataset avoids leakage caused by highly correlated adjacent depth samples. The training set consists of five wells, while separate wells are used for validation and testing (Table 2). This partitioning reflects the practical deployment scenario, where models are required to generalize to unseen wells with different lithological distributions.

A representative well profile is shown in Figure 1, illustrating the correspondence between logging responses and lithological variations. This visualization highlights the core challenge addressed in this study: similar logging signatures may correspond to different mechanical properties depending on lithological regime, confirming the non-uniqueness of the underlying mapping.

In summary, the constructed dataset is specifically designed to expose and preserve lithotype-induced conditional variability. This enables the proposed model to learn structured residuals associated with heterogeneous geological conditions, rather than fitting a single global mapping. Consequently, the dataset provides a rigorous foundation for evaluating whether the proposed method can help address heterogeneity-induced non-uniqueness and restore physically consistent geomechanical characterization.

2.3. Lithotype Classification and Feature Engineering Based on the HMLZ Index

2.3.1. Experimental Determination of Rock Strength Parameters

To establish a reliable reference for evaluating the proposed heterogeneity-aware modeling framework, laboratory geomechanical experiments were conducted on coal and parting (dirt band) samples collected from deep coal seams in the Ordos Basin. The objective of these experiments is not only to provide ground-truth labels for supervised learning, but more importantly to verify whether the proposed model can recover physically consistent geomechanical responses under heterogeneous lithological conditions.

Key mechanical parameters, including Uniaxial Compressive Strength (UCS), Elastic Modulus (E), and Poisson’s ratio ($\nu$), were measured. These measurements constitute an independent experimental benchmark, enabling direct comparison between predicted and observed values at corresponding depths. This design ensures that model evaluation is grounded in physical measurements rather than solely relying on statistical error metrics, thereby providing a rigorous basis for assessing whether the model resolves heterogeneity-induced bias.

The experiments were performed using the MTS-816 electro-hydraulic servo-controlled rock testing system (Figure 2a), which provides stable closed-loop control for load application and displacement measurement. Uniaxial compression tests were conducted under zero confining pressure until specimen failure. The UCS was calculated from the peak load:

$$UCS={\sigma }_c={\displaystyle \frac{F_{max}}{A}}$$

(2)

where A is the initial cross-sectional area. The elastic modulus E was determined from the slope of the stress–strain curve within the linear elastic regime, while Poisson’s ratio $\nu$ was computed as the ratio of lateral strain to axial strain in the same region.

During testing, failure modes were recorded to capture lithotype-dependent structural behavior. These observations provide qualitative evidence of mechanical variability across lithological regimes and support the interpretation that geomechanical response is conditioned by lithotype rather than solely by logging signatures.

Specimen preparation followed standard rock mechanics protocols to ensure comparability. Cylindrical samples were fabricated with strict control of geometry, and end faces were precision-ground to minimize loading artifacts. Defective specimens were excluded to maintain data integrity. Special care was taken during handling and mounting to avoid artificial damage, particularly given the fragile and heterogeneous nature of coal.

Figure 3. Measured UCS, elastic modulus, and Poisson’s ratio for samples in Well-Train-3.

Figure 3. Measured UCS, elastic modulus, and Poisson’s ratio for samples in Well-Train-3.

Experimental results reveal clear lithotype-dependent patterns. Coal samples exhibit relatively low strength and stiffness, with UCS values ranging from 7.7 to 20.5 MPa and elastic modulus between 4.1 and 7.8 GPa, while partings show significantly higher strength (up to 71.8 MPa). Poisson’s ratio is generally higher in coal (0.23–0.34), reflecting its deformable structure. Moreover, substantial variability exists among different coal lithotypes, confirming that mechanical response is strongly conditioned by lithological regime.

These observations directly support the central premise of this study: identical or similar logging responses may correspond to different mechanical properties depending on lithotype, leading to a non-unique mapping. Therefore, lithotype-induced variability should not be treated as noise, but as a structured component of the geomechanical response.

The measured parameters are subsequently used as reference targets for model evaluation. Predictions at corresponding depths are extracted and compared against experimental values to assess both accuracy and physical consistency under cross-well conditions. This validation strategy enables a direct examination of whether the proposed model successfully captures lithotype-conditioned deviations rather than merely fitting global trends.

In summary, the experimental dataset provides both quantitative supervision and physical evidence of heterogeneity-induced non-uniqueness. This forms a critical foundation for evaluating the proposed structured learning formulation, in which geomechanical response is decomposed into global and lithotype-dependent components.

Figure 4. Measured mechanical properties for samples in Well-Train-5.

Figure 4. Measured mechanical properties for samples in Well-Train-5.

2.3.2. Feature Construction

Based on the experimentally measured rock mechanical parameters, the input representation is constructed to support the proposed heterogeneity-aware residual learning paradigm. The objective is not merely to assemble predictive features, but to establish a structured input space that enables the modeling of lithotype-conditioned deviations in geomechanical response.

The input is defined as a triplet $(X,L,z)$, where X denotes logging-derived physical measurements, L represents lithotype as a lithotype condition, and z corresponds to depth. This formulation reflects the assumption that geomechanical response is governed by both global trends encoded in X and structured, lithotype-dependent variations captured through L.

The primary feature set X consists of conventional well logging curves, including acoustic transit time (AC), bulk density (DEN), gamma ray (GR), compensated neutron logging (CNL), and resistivity measurements (LLD, LLS). These features collectively describe formation properties from multiple physical perspectives, including elastic response, density distribution, mineral composition, porosity, and fluid-related conductivity. Importantly, these measurements alone do not uniquely determine mechanical properties under heterogeneous conditions, as similar logging signatures may correspond to different lithological regimes.

Depth (z) is incorporated as a continuous variable to preserve vertical continuity and capture systematic variations associated with burial conditions and stress environments. This inclusion ensures that large-scale geological trends are retained in the global component of the model.

To explicitly account for lithological heterogeneity, lithotype information derived from the HMLZ index is introduced as a lithotype condition L. Unlike conventional feature-augmented approaches, lithotype is not treated as an additional predictor within a single mapping. Instead, it governs the residual component of the model, enabling the learning of lithotype-dependent corrections to the global response. This design directly addresses cases where identical logging responses correspond to different mechanical properties, thereby helping address heterogeneity-induced non-uniqueness.

Figure 5. Pearson correlation coefficient matrix between logging features and rock mechanical parameters.

Figure 5. Pearson correlation coefficient matrix between logging features and rock mechanical parameters.

To assess the relevance of the selected features, Pearson correlation coefficients between logging variables and mechanical parameters are computed [16]. The results show that the selected logging features exhibit meaningful correlations with the target variables, indicating that they provide informative signals for geomechanical characterization. However, these correlations are primarily linear and do not fully capture the complex, nonlinear, and lithotype-dependent relationships present in the data.

Therefore, while the feature set provides a physically meaningful basis for modeling, the key improvement of the proposed method does not stem from feature selection itself. Instead, it arises from the structured modeling of conditional bias through lithotype-conditioned residual learning. In this context, the input representation $(X,L,z)$ serves as the foundation for decomposing geomechanical response into global and lithotype-dependent components.

2.3.3. Coal Lithotype Identification (HMLZ)

Coal seams exhibit strong lithological heterogeneity governed by depositional environment, maceral composition, and structural development. Variations in fracture density, pore structure, and compositional characteristics lead to substantial differences in geomechanical behavior, even under similar logging responses [17]. This heterogeneity induces non-uniqueness in the mapping between logging data and mechanical properties, which cannot be resolved by models assuming a single global relationship.

To explicitly represent lithological regimes and enable conditional modeling of geomechanical response, the HMLZ index is employed to construct a lithotype-dependent representation along the wellbore. Rather than serving merely as a classification tool, the HMLZ-derived lithotype sequence is used to represent variations in mechanical behavior under heterogeneous conditions.

The HMLZ index is defined based on conventional well logging measurements, including resistivity (RD), acoustic transit time (AC), bulk density (DEN), and gamma ray (GR), capturing key petrophysical characteristics associated with coal brittleness:

$$HMLZ={\displaystyle \frac{lg\left(RD\right)\times AC}{{DEN}^2\times GR}}$$

(3)

Based on the computed HMLZ values, coal lithotypes are categorized into four classes—bright coal, semi-bright coal, semi-dull coal, and dull coal—according to predefined threshold intervals. These lithotypes reflect systematic differences in maceral composition, structural integrity, and fracture development, which are directly linked to mechanical properties. The classification criteria and corresponding characteristics are summarized in Table 3.

The HMLZ index is computed continuously along depth, producing a lithotype sequence $L\left(z\right)$ that characterizes the spatial distribution of heterogeneity. Within the proposed framework, this sequence is treated as a lithotype condition rather than an input feature in a unified mapping. Specifically, $L\left(z\right)$ governs the residual component of the model, enabling the learning of lithotype-dependent deviations from the global geomechanical response.

This design directly addresses the key challenge of heterogeneous systems: identical logging responses may correspond to different mechanical properties under different lithological regimes. By conditioning the residual model on lithotype, the framework captures structured, regime-dependent variations that would otherwise be averaged out in conventional approaches.

Importantly, the performance gain achieved by incorporating HMLZ does not arise from additional information alone, but from the structural modeling of conditional bias. In this sense, the HMLZ-derived lithotype sequence provides a physically meaningful partition of the data space, enabling the decomposition of geomechanical response into global and lithotype-dependent components.

In summary, the HMLZ-based lithotype identification establishes a critical link between geological heterogeneity and model structure. It enables the proposed method to move toward a structured learning formulation that explicitly accounts for lithotype-controlled variations in geomechanical behavior.

2.4. Lithotype-Conditioned Residual Characterization Framework

To improve geomechanical characterization under lithological heterogeneity, a lithotype-conditioned residual learning framework is established in this study. The underlying motivation is that, in heterogeneous coal-bearing formations, similar logging responses may correspond to different mechanical properties under different lithological regimes. Under such conditions, a single global mapping is often insufficient to fully describe the response relationship. To address this issue, the prediction is formulated as the sum of a global component and a lithotype-conditioned residual component.

Accordingly, the final prediction is expressed as

$$\widehat{Y}\left(d\right)=f\left(X\left(d\right)\right)+g(X\left(d\right),L\left(d\right))$$

(4)

where $X\left(d\right)$ denotes the logging-derived input features at depth d, $L\left(d\right)$ denotes the lithotype label derived from the HMLZ index, $f(·)$ represents the global mapping from logging responses to mechanical parameters, and $g(·)$ represents the lithotype-conditioned residual correction.

In implementation, both $f(·)$ and $g(·)$ are constructed using CatBoost regressors. The global model $f\left(X\right)$ is trained using only logging features, while the residual model $g(X,L)$ takes both logging features and lithotype as inputs. Lithotype is encoded as a categorical variable within CatBoost.

For multi-output prediction, separate models are trained for each target variable to ensure stable optimization. Hyperparameters are optimized using Bayesian Optimization based on validation error, and the same optimization protocol is applied to both stages to ensure fair comparison.

Residuals used to train $g(·)$ are computed on the training set using predictions from $f(·)$ without data leakage.

The function $f\left(X\right)$ is used to learn the dominant relationship between logging responses and geomechanical parameters over the entire training dataset. It captures the general response trend shared across samples and reflects the lithology-independent component of the prediction. However, because the data are affected by lithological heterogeneity, this global model alone may leave systematic errors in intervals where different lithotypes exhibit different mechanical responses under similar logging signatures.

To characterize this effect explicitly, the residual is defined as

$$r\left(d\right)=Y\left(d\right)-f\left(X\right(d\left)\right)$$

(5)

where $Y\left(d\right)$ is the measured target vector and $r\left(d\right)$ is the deviation between the observation and the global prediction. In the proposed framework, this residual is interpreted as a structured correction term associated with lithological heterogeneity rather than as purely random noise. The function $g(X,L)$ is then used to learn the relationship between the residual, the logging responses, and the lithotype condition.

The training procedure is implemented sequentially. First, the global model $f\left(X\right)$ is trained using the logging features to predict the target mechanical parameters. Second, the residuals are computed on the training set as the difference between the measured values and the predictions of the global model. Third, a residual model $g(X,L)$ is trained using the logging features together with the lithotype variable, with the residual term as the prediction target. During inference, the output of the global model and that of the residual model are added to obtain the final prediction.

In practical implementation, lithotype is introduced as a categorical lithotype condition in the residual stage rather than being used only as an ordinary feature in a single unified predictor. This design allows the residual model to learn lithotype-dependent corrections to the global trend and thereby improves adaptability in heterogeneous intervals. The formulation does not assume that all samples follow exactly the same response relationship; instead, it allows systematic deviations associated with lithological regime to be represented in an explicit manner.

This decomposition also improves the interpretability of the modeling framework. The global component describes the dominant mapping shared by the dataset, whereas the residual component accounts for lithotype-related deviations from that common trend. In this sense, the final prediction can be understood as a combination of baseline geomechanical response and lithology-dependent correction.

It should be noted that the present framework is developed and validated using data from coal-bearing formations in the Ordos Basin. The applicability of this approach to other basins, lithologies, and logging configurations requires further verification.

In summary, the proposed lithotype-conditioned residual framework reformulates geomechanical characterization under heterogeneity as a decomposition problem composed of a global predictor and a lithotype-dependent correction term. This provides a more explicit way to represent heterogeneity-induced deviations and offers a physically more consistent basis for prediction under the geological conditions considered in this study.

2.4.1. Hyperparameter Optimization and SHAP Analysis

To ensure stable model construction, hyperparameter optimization was performed using Bayesian Optimization (BO). In this study, BO was used to search for suitable parameter combinations for the predictive model by minimizing the Mean Squared Error on the validation set. The optimized parameters mainly include the number of iterations, learning rate, tree depth, and L2 regularization coefficient. This procedure improves the stability of model training and reduces the risk of overfitting caused by manual parameter selection.

The purpose of hyperparameter optimization in this study is to obtain a stable and reproducible model configuration for subsequent comparison and analysis. It should be emphasized that the performance improvement of the proposed framework is not attributed to hyperparameter tuning itself, but to the introduction of lithotype-conditioned residual modeling under the same optimization protocol.

In addition, SHAP analysis was employed to examine the contribution patterns of the input features in the trained model. SHAP provides an additive explanation of model prediction and is used here only for supplementary interpretation rather than as primary evidence of the proposed formulation. It enables the contribution of each feature to be quantified at both the global and sample levels. For a given sample, the SHAP formulation can be written as

$$g\left(z^{\prime }\right)={\varphi }_0+\sum _{j=1}^M{\varphi }_j{z}_j^{\prime }$$

(6)

where ${\varphi }_0$ is the baseline output and ${\varphi }_j$ represents the contribution of the jth feature.

The SHAP value of each feature is computed as the weighted average marginal contribution over all possible feature subsets:

$${\psi }_j=\sum _{S\subseteq \{x_1,\dots ,x_p\}\setminus \left\{x_j\right\}}{\displaystyle \frac{\left|S\right|!(p-|S|-1)!}{p!}}\left(f_x(S\cup \left\{x_j\right\})-f_x\left(S\right)\right)$$

(7)

In this work, TreeSHAP was used for efficient explanation of the tree-based model. The SHAP analysis was mainly used as an auxiliary interpretive tool to examine whether the learned prediction behavior is consistent with known geomechanical understanding. Specifically, global SHAP importance was used to identify the dominant logging variables affecting the prediction, while local SHAP analysis was used to inspect how feature contributions vary in different depth intervals and lithological settings.

It should also be noted that SHAP is not used here as the primary evidence for validating the proposed residual formulation. The main evidence for the effectiveness of the framework is still derived from cross-well evaluation, case-study comparison, and ablation analysis. SHAP is used only to provide supplementary interpretation of the trained model and to help assess whether the learned feature-response relationships remain physically plausible.

Therefore, BO and SHAP play different roles in the present study. BO is used to improve the stability of model training, whereas SHAP is used to provide auxiliary interpretive support for model behavior. Together, they complement the quantitative evaluation of the proposed lithotype-conditioned residual framework.

3. Results and Discussion

3.1. Reliability-Oriented Hyperparameter Optimization and Model Performance

To obtain a stable training configuration for the proposed model under heterogeneous data conditions, hyperparameter optimization was conducted using Bayesian Optimization (BO) in combination with cross-validation (CV) [18]. The objective of this procedure is to identify a reliable set of hyperparameters that minimizes validation error while maintaining stable model behavior across different data splits.

BO iteratively refines the hyperparameter configuration by modeling the validation error as a black-box function and selecting candidate parameter sets through an acquisition strategy. Compared with manual tuning or grid-based search, this approach enables more efficient exploration of the parameter space and provides a more reproducible model-selection process.

As shown in Figure 6, the validation error exhibits clear trends with respect to the tested hyperparameters. Increasing the number of iterations reduces the validation error up to a certain point, after which the improvement becomes marginal, indicating convergence in model training. The learning rate also shows a limited effective range: excessively small values lead to underfitting, whereas overly large values result in unstable training and poorer validation performance.

The optimal tree depth is relatively low, suggesting that a simple model structure is sufficient to capture the dominant relationships in the data. Increasing the depth leads to higher validation error, indicating overfitting and increased sensitivity to noise. Similarly, larger values of the regularization coefficient (l2_leaf_reg) tend to increase the validation error, implying that excessive regularization may weaken the model’s ability to fit the data adequately.

These observations indicate that the prediction task favors a configuration with moderate model complexity and stable training behavior. Overall, the selected hyperparameter region is consistent with the need to balance fitting capacity and regularization in a heterogeneous prediction setting.

Table 4. Optimal CatBoost hyperparameters obtained after tuning.

Parameter	Description	Search Range	Optimal Value
iterations	Number of boosting iterations	[40, 200]	150
learning rate	Step size controlling update magnitude	[0.01, 0.5]	0.38
depth	Maximum depth of decision trees	[2, 10]	3
l2_leaf_reg	L2 regularization coefficient	[0.01, 1]	0.04

Table 4. Optimal CatBoost hyperparameters obtained after tuning.

Parameter	Description	Search Range	Optimal Value
iterations	Number of boosting iterations	[40, 200]	150
learning rate	Step size controlling update magnitude	[0.01, 0.5]	0.38
depth	Maximum depth of decision trees	[2, 10]	3
l2_leaf_reg	L2 regularization coefficient	[0.01, 1]	0.04

Based on the optimization results, the final hyperparameter configuration is set as iterations = 150, learning rate = 0.38, depth = 3, and l2_leaf_reg = 0.04. This configuration provides a practical balance between model capacity and regularization, leading to stable training behavior in subsequent experiments.

Importantly, hyperparameter optimization should be viewed as a supporting training procedure rather than the primary source of performance improvement. Its role is to provide a stable and consistent experimental configuration, whereas the comparative performance gains are evaluated in the following sections with respect to the proposed modeling formulation.

In summary, hyperparameter optimization serves as a supporting training procedure that helps establish a stable and reproducible model configuration, rather than being the primary source of performance gain.

3.2. Lithotype-Aware Characterization of Rock Mechanical Parameters

Following data preprocessing and model optimization, the proposed heterogeneity-aware residual formulation was applied to the prediction wells for continuous characterization of rock mechanical parameters. Table 5 summarizes the performance across the training, validation, and test sets. The model achieves an $R^2$ of 0.982 and a Mean Absolute Error (MAE) of 0.312 MPa on the training set, indicating strong fitting capacity. On the validation and test sets, the $R^2$ values are 0.936 and 0.928, with MAE values of 0.594 and 0.641 MPa, respectively.

The consistency of these metrics across datasets indicates that the model maintains stable predictive behavior under cross-well conditions, without evident overfitting.

To further examine the effectiveness of the method under heterogeneous conditions, Well-Validate-1 and Well-Test-1 were selected as representative cases for comparison with conventional empirical approaches. Traditional methods assume lithological homogeneity and apply a single parametric relationship across the entire well interval. Under heterogeneous conditions, this assumption leads to systematic errors, particularly in lithological transition zones where identical logging responses may correspond to different mechanical properties.

In contrast, the proposed formulation separates geomechanical response into a global component and lithotype-dependent deviations. This enables the model to adapt its prediction according to lithological regime, capturing variations at lithological interfaces and reducing systematic bias in heterogeneous intervals.

The observed performance improvement should not be interpreted solely as an increase in prediction accuracy. Instead, it reflects a mitigation of the non-uniqueness in the mapping between logging responses and mechanical properties under heterogeneous conditions. By explicitly modeling lithotype-conditioned residuals, the method provides additional flexibility to distinguish cases that are difficult to separate under a single global mapping.

Furthermore, the improvement is not solely attributable to the inclusion of additional information. The same logging features are used in both conventional and proposed approaches; the difference lies in how the mapping is structured. The proposed method explicitly models lithology-related variations, whereas conventional methods implicitly average over such effects.

As a result, the method reduces lithotype-induced bias and improves the consistency between predicted mechanical parameters and underlying geological conditions. This provides a more reliable basis for geomechanical characterization in heterogeneous formations and suggests that the performance gain is associated with the proposed modeling formulation.

In summary, the results indicate that the proposed approach helps mitigate the non-uniqueness in geomechanical characterization under heterogeneous conditions, while maintaining stable predictive performance and physical consistency.

3.2.1. Case Study: Well-Validate-1

Well-Validate-1 was selected as a representative validation well to examine the behavior of the proposed method under heterogeneous lithological conditions. In the depth interval of 2700–2950 m, a relatively stable coal seam with a thickness of approximately 48 m is developed, interbedded with 4–6 parting layers totaling about 12 m. A total of 23 experimental measurements are available within this interval, covering all four lithotypes (bright, semi-bright, semi-dull, and dull coal), providing a suitable basis for evaluating both predictive accuracy and lithotype-dependent behavior. The predicted profiles, including HMLZ, lithotype classification, and mechanical parameters, are shown in Figure 7.

The predicted UCS profile shows clear alignment with lithotype variations along the depth axis. In the bright coal interval (2750–2800 m), where HMLZ values indicate low-strength lithotypes, the predicted UCS remains within a low range (18–28 MPa), consistent with measured values and yielding a Mean Absolute Error (MAE) of 2.3 MPa. Local fluctuations within this interval are also captured, reflecting sensitivity to small-scale structural variations.

At the lithological transition near 2800 m, an increase in HMLZ corresponds to a corresponding rise in predicted UCS over a short depth interval. This transition is closely aligned with measured data, indicating that the model responds to lithotype changes rather than producing a smoothed global trend. In the dull coal section around 2820 m, the predicted UCS stabilizes within the high-strength range and remains consistent with experimental observations.

In contrast, the traditional empirical approach exhibits systematic distortion across lithotypes. In bright coal, it overestimates UCS due to bias toward higher-strength samples, while in dull coal it underestimates strength. More critically, the traditional prediction curve lacks sensitivity to lithological transitions, resulting in a smooth, low-frequency trend that fails to capture abrupt changes in mechanical behavior.

A lithotype-wise comparison reveals that errors in the traditional method are concentrated in extreme lithotypes, whereas the proposed method maintains consistently low errors across all categories. Across all 23 samples, the proposed method reduces the MAE from 9.2 MPa to 3.1 MPa and improves the coefficient of determination from 0.69 to 0.94.

These results should not be interpreted solely as an improvement in prediction accuracy. Instead, they indicate that the proposed formulation helps reduce lithotype-induced bias by explicitly modeling conditional deviations. The key distinction is that the method does not enforce a single global mapping, but adapts predictions according to lithological regime.

Further evidence of this mechanism can be observed in the behavior of other mechanical parameters. The elastic modulus decreases in low-strength lithotypes and increases in high-strength lithotypes, while Poisson’s ratio exhibits the opposite trend. These patterns are consistently captured by the proposed method but are systematically distorted in the traditional approach. The agreement across multiple parameters indicates that the model preserves physically consistent relationships rather than fitting isolated targets.

Overall, the case study demonstrates that lithotype-dependent deviations are not random noise but structured variations governed by geological conditions. By explicitly modeling these variations, the proposed method captures regime-dependent behavior and helps mitigate the non-uniqueness inherent in heterogeneous systems. This suggests that the performance gain is associated with the modeling of lithology-dependent variations.

3.2.2. Case Study: Well-Test-1

Well-Test-1 serves as an independent test well that was not involved in either the training or validation stages. It is therefore used to evaluate the generalization capability of the proposed heterogeneity-aware formulation under more complex geological conditions. Compared with Well-Validate-1, this well exhibits significantly stronger lithological heterogeneity. Within the 2760–2920 m interval, the HMLZ curve identifies 37 lithotype transitions, corresponding to an average spacing of approximately 4.3 m, indicating a high-frequency heterogeneous system. In addition, multiple thin parting layers are present in the 2800–2850 m interval, forming frequent interbedded contacts with coal seams. This configuration imposes stringent requirements on the model’s ability to resolve rapid lithological transitions and to maintain consistency across narrow depth intervals. A total of 18 experimental measurements are available, including nine points located within high-frequency transition zones, providing a demanding test scenario. The prediction results are shown in Figure 8.

In the high-frequency transition interval (2800–2850 m), the proposed method shows clear responsiveness to lithotype variations. Multiple step-like changes in UCS are captured within short depth ranges, and the prediction curve remains synchronized with lithotype transitions indicated by the HMLZ profile. At lithological interfaces, the model produces rapid and localized adjustments in predicted strength, consistent with measured data. Thin parting layers are also correctly identified as high-strength zones, with predictions transitioning smoothly to adjacent coal intervals.

In contrast, the traditional empirical approach shows a loss of stability under these conditions. The prediction curve exhibits irregular oscillations that are not aligned with lithotype changes, and abrupt variations appear even in relatively stable intervals. This behavior is consistent with the limitation of enforcing a single global mapping, which cannot accommodate rapid regime changes.

Quantitative evaluation further highlights this difference. In the high-frequency transition zone, the proposed method maintains low prediction error and a high proportion of samples within a narrow error band, whereas the traditional method shows significantly larger deviations and reduced consistency. Across the entire well, the proposed method achieves substantially lower error and higher correlation with measured values, while maintaining similar performance levels to those observed in the validation well. The small difference in error between validation and test wells suggests that the learned mapping maintains stable predictive performance on unseen data.

These results should not be interpreted solely as accuracy improvement. Under high-frequency heterogeneity, the mapping from logging responses to mechanical properties becomes highly non-unique, and the residual component becomes more significant. The ability of the proposed method to maintain stable predictions suggests that the residual term captures lithotype-dependent variations rather than random fluctuations.

In particular, the alignment between predicted step changes and lithotype transitions suggests that the residual is not purely noise, but is related to lithological regime. This indicates that the decomposition into global and lithotype-conditioned components provides a reasonable representation of the underlying geomechanical behavior.

Overall, the comparison across validation and test wells shows that the proposed formulation remains stable under both moderate and high-frequency heterogeneity. The improvement in performance is therefore associated with the explicit modeling of lithological-regime-related variations, rather than with increased model complexity or additional features.

This suggests that the proposed approach helps mitigate the non-uniqueness in geomechanical characterization under heterogeneous conditions, achieving consistent cross-well generalization while preserving physical interpretability.

3.3. Ablation Study Results and Quantitative Analysis

To quantify the role of the lithotype-conditioned residual formulation, three modeling configurations were evaluated on the independent test well (Well-Test-1). The comparison focuses on uniaxial compressive strength (UCS) prediction and is summarized in Table 6.

The baseline model, which relies solely on logging features, captures the overall variation trend but exhibits significant errors in intervals where similar logging responses correspond to different lithotypes. This suggests that a single global mapping may be insufficient under heterogeneous conditions.

Introducing lithotype as an additional input improves performance, indicating that lithological information provides useful constraints. However, the improvement remains limited, suggesting that treating lithotype as a conventional feature does not fully resolve the ambiguity in the mapping.

In contrast, the proposed formulation models lithotype-induced deviations as a separate residual component. This leads to a noticeable performance improvement, with $R^2$ increasing to 0.928 and MAE reduced by more than 60% compared to the baseline. More importantly, this improvement is achieved without introducing new information, but by restructuring the mapping itself.

The comparison between the feature-augmented model and the proposed formulation provides evidence that the residual is not purely random noise. If lithotype-dependent variations were unstructured, incorporating lithotype as a feature would be sufficient. However, the additional improvement achieved by the residual formulation indicates that these variations may follow a systematic pattern that cannot be captured within a single mapping.

This result suggests that lithotype is associated with condition-dependent deviations in geomechanical response. By explicitly modeling these deviations, the proposed method separates global trends from lithotype-specific corrections, thereby helping mitigate the non-uniqueness inherent in heterogeneous systems.

In this sense, the performance gain should be interpreted as a consequence of modeling conditional bias rather than improving predictive capacity. The results suggest that the relationship between logging responses and mechanical parameters may exhibit multi-regime characteristics, and that accurate characterization requires explicit representation of lithotype-dependent structure.

Overall, the ablation study provides quantitative evidence that the proposed approach differs from conventional formulations and is a structurally different formulation that addresses limitations associated with global mapping assumptions in heterogeneous geological environments.

3.4. Mechanism Analysis of Structured Residual Correction

To verify that the performance gain of the proposed framework arises from modeling structured lithotype-dependent deviations rather than from a purely empirical two-stage refinement, a dedicated residual analysis was conducted. This analysis examines the error distribution of the baseline global model $f\left(X\right)$ and evaluates how the residual component $g(X,L)$ systematically resolves heterogeneity-induced bias.

3.4.1. Identification of Structured Bias in the Global Baseline

A fundamental assumption of this study is that a single global mapping $f\left(X\right)$ inevitably introduces systematic bias in heterogeneous formations because it “averages” the mechanical responses of different lithotypes. To test this, the Signed Mean Residual (SMR) and Standard Deviation (SD) of the baseline model were calculated for each coal lithotype in the test set.

As shown in Table 7, the baseline residuals are not randomly distributed white noise; instead, they exhibit a clear polarity tied to the lithological regime. In low-strength bright coal intervals, the baseline model consistently overestimates UCS (SMR = +2.45 MPa), whereas in high-strength dull coal intervals, it tends to underestimate the strength (SMR = -3.12 MPa). This systematic departure suggests that the “ambiguity” mentioned in the introduction—where similar logging responses correspond to different mechanical properties—manifests as a structured bias in a unified predictor.

3.4.2. Amplification of Bias in Transition Zones

The structured bias is further intensified in lithological transition zones. We defined “Transition Zones” as intervals within 0.5 m of a lithotype boundary identified by the HMLZ index. Figure 9 compares the residual density between stable intervals and transition zones.

In stable lithological intervals, the baseline model exhibits a relatively narrow error distribution. However, in transition zones, the residual variance increases by approximately 140%, and the distribution becomes markedly bimodal. This phenomenon indicates that near lithological interfaces, the global mapping fails to track the rapid shift in geomechanical response, even when the logging signals (e.g., AC or DEN) show only subtle variations. This is consistent with the conclusion that heterogeneity-induced ambiguity is a localized stress-point for traditional modeling approaches.

3.4.3. Effectiveness of the Lithotype-Conditioned Correction

The proposed method addresses this by explicitly modeling these structured errors through $g(X,L)$. Figure 10 illustrates the “flattening” effect of the residual correction. After incorporating the lithotype-conditioned component, the SMR for all lithotypes converged toward zero (e.g., dull coal SMR improved from -3.12 MPa to -0.18 MPa).

Crucially, the standard deviation of the residuals also decreased across all regimes, indicating that $g(X,L)$ does not just shift the mean but also reduces the uncertainty within each lithotype. This transition from a “lithotype-biased” error to a “lithotype-neutral” error provides strong evidence that the performance gain is associated with the proposed formulation. The residual model successfully captures the conditional deviations induced by the HMLZ-defined regimes, improving the physical consistency of the characterization.

This observation indicates that the residual is not only structured, but also explicitly dependent on lithotype, confirming that lithological regimes act as conditioning variables governing systematic deviations in geomechanical response. This supports the assumption that the mapping from X to Y is multi-regime rather than globally unique.

In summary, the residual analysis demonstrates that: (1) baseline errors are geologically structured; (2) this structure is driven by lithological heterogeneity; and (3) the proposed decomposition captures this structure, transforming a biased global mapping into a lithotype-aware characterization framework.

3.5. Heterogeneity-Focused Evaluation in Transition Zones

To further examine the role of lithotype-induced heterogeneity, a focused evaluation was conducted on the 2800–2850 m interval of Well-Test-1, where bright coal and dull coal are frequently interbedded. This interval represents a typical high-frequency heterogeneous regime, in which the mapping between logging responses and mechanical properties becomes highly non-unique. The quantitative results are summarized in Table 8.

The baseline model shows a pronounced degradation in performance within this interval, producing overly smoothed predictions that fail to capture sharp variations in mechanical properties across lithological boundaries. This behavior is consistent with the limitation of a single global mapping when applied to a multi-regime system.

Introducing lithotype as an additional feature improves sensitivity to lithological variation, but substantial errors remain. This indicates that while lithotype contains relevant information, embedding it within a unified mapping does not sufficiently resolve the ambiguity caused by heterogeneity.

In contrast, the proposed formulation maintains stable predictive performance in the transition zone. The improvement is particularly evident in the reduction of maximum error, suggesting that abrupt changes in mechanical properties are better captured. This behavior suggests that the model is able to adapt its response locally in accordance with lithological transitions.

More importantly, this interval provides additional evidence of the role of the residual component. In transition zones, the global mapping becomes insufficient and the residual term becomes more significant. The significant performance gap between the feature-augmented model and the proposed formulation indicates that lithotype-induced deviations are not purely random and require explicit modeling rather than implicit representation.

This observation supports the interpretation that heterogeneity introduces condition-dependent bias into the mapping. By modeling this bias as a lithotype-conditioned residual, the proposed method helps mitigate the ambiguity that arises in transition zones and maintains consistent predictive behavior.

Overall, the results demonstrate that the advantage of the proposed formulation is most pronounced in intervals where heterogeneity is strongest. This suggests that the method addresses limitations of global mapping approaches and provides a reliable characterization of mechanical properties in complex geological settings.

In this sense, transition zones serve as a critical test for evaluating whether the residual component captures structured variations. The consistent improvement observed in this interval indicates that the residual is not noise, but is related to lithological regime, providing support for the proposed heterogeneity-aware modeling framework.

3.6. SHAP-Based Interpretability Analysis

To examine whether the proposed heterogeneity-aware formulation captures physically meaningful and lithotype-consistent behavior under heterogeneous geological conditions, SHAP (SHapley Additive Explanations) was employed to analyze the relationships between logging responses, lithological regimes, and predicted mechanical parameters. In this study, SHAP is not used merely as a post hoc explanation tool, but as an analysis tool to examine model behavior.

At the global level, the mean absolute SHAP value was calculated for both the training and test sets to rank feature importance. As shown in Table 9, acoustic transit time (AC) is the most influential feature in both datasets, followed by gamma ray (GR), density (DEN), and resistivity-related features. The overall consistency of the feature ranking between the training and test sets indicates that the model captures stable controlling factors across wells, rather than overfitting to specific local patterns.

This global ranking is physically interpretable. AC reflects the propagation behavior of acoustic waves and is closely associated with fracture development, pore structure, and structural integrity. DEN characterizes material compactness and bulk structural condition, while GR provides supplementary information related to compositional variability, ash content, and clay-related effects. Together, these features form a physically meaningful basis for geomechanical characterization.

To further examine the direction and distribution of feature effects, a SHAP summary plot was generated, as shown in Figure 11. The summary plot presents the SHAP contribution of each sample together with the corresponding feature value, thereby revealing both the magnitude and polarity of feature influence. High AC values generally correspond to negative SHAP contributions, indicating a reduction in the predicted mechanical parameters, whereas low AC values tend to produce positive contributions. This agrees with rock physics expectations, since larger transit time is usually associated with poorer structural integrity and lower load-bearing capacity. In contrast, higher DEN values generally produce positive contributions, reflecting the higher strength expected in denser and more compact formations.

Unlike simple correlation analysis, SHAP can reveal conditional effects arising from nonlinear interactions and lithotype-dependent modulation. To assess whether lithological heterogeneity is explicitly reflected in the model behavior, the SHAP contributions associated with coal lithotypes were further examined. The analysis shows that even within similar AC or DEN intervals, SHAP values exhibit systematic offsets across lithotypes. Lower-strength lithotypes tend to produce more negative contributions, whereas higher-strength lithotypes more often produce positive contributions. This indicates that identical logging responses do not correspond to a unique mechanical implication, but are interpreted differently depending on lithological regime.

This result provides additional support for the proposed formulation. If lithotype-induced variation were merely random noise, samples with similar logging features would show similar SHAP contributions regardless of lithotype. Instead, the observed systematic separation suggests that the model captures regime-dependent behavior, which is consistent with the assumption that heterogeneity introduces condition-dependent variations into the mapping from logging data to mechanical properties.

To visualize this modulation more directly, local contribution analyses were performed for representative depth samples. Figure 12 shows how the final prediction is decomposed into additive feature contributions for specific examples. For low-strength samples located in bright coal intervals, the prediction below the baseline is typically driven by the joint effect of high AC, low DEN, and lithotype-associated negative contributions. For high-strength samples in dull or semi-dull coal intervals, the opposite pattern is observed, with low AC, high DEN, and lithotype-associated positive contributions driving the prediction above the baseline. This decomposition provides a traceable explanation for why the predicted mechanical parameter at a given depth is high or low.

Taken together, the SHAP results provide consistent evidence from multiple perspectives: global importance ranking, direction of feature influence, lithotype-dependent modulation, and local sample-wise decomposition. More importantly, they support the central claim of this study that the performance gain does not arise from feature enrichment alone, but from the structural modeling of lithotype-conditioned deviations. The interpretability analysis therefore suggests that the proposed method captures physically meaningful and geologically consistent patterns, and that the residual component represents lithotype-dependent behavior rather than arbitrary correction.

In summary, the SHAP-based analysis demonstrates that the proposed method preserves physical consistency, reflects lithotype-controlled modulation of geomechanical response, and provides additional interpretability support for the heterogeneity-aware residual formulation. These findings reinforce the conclusion that the model helps mitigate heterogeneity-induced non-uniqueness through structured conditional modeling rather than through incremental adjustment of a single global predictor.

4. Conclusions

This study addresses the challenge of geomechanical characterization under lithological heterogeneity, where similar logging responses may correspond to different mechanical properties, leading to ambiguity in the mapping between input and target variables. The results indicate that performance improvement is primarily associated with explicitly modeling heterogeneity-induced deviations, rather than increasing model complexity or incorporating additional features.

A lithotype-aware residual learning formulation is introduced to decompose geomechanical response into a global component and lithotype-dependent deviations. This formulation enables the representation of conditional bias associated with heterogeneous lithological regimes and provides a structured alternative to conventional single-mapping approaches.

Compared with traditional methods that assume a unified mapping between logging responses and mechanical parameters, the proposed approach captures lithotype-dependent variations through residual modeling. This reduces systematic bias in heterogeneous intervals and improves the consistency between predicted mechanical properties and underlying geological conditions.

Cross-well evaluation shows that the method maintains stable performance across training, validation, and independent test wells, suggesting that the learned relationships are not restricted to individual wells. Case studies further demonstrate that the model responds effectively to lithological transitions and reduces prediction errors in extreme lithotypes, particularly in intervals where heterogeneity is pronounced.

Ablation experiments indicate that incorporating lithotype as a standard feature provides limited improvement, whereas residual decomposition introduces additional gains by modeling structured deviations. Interpretability analysis suggests that feature contributions exhibit lithotype-dependent variation, supporting the role of the residual component in capturing conditional behavior.

Despite these findings, several limitations should be noted. The dataset is derived from a specific coal-bearing formation within a single basin, and the number of wells is limited. Therefore, the applicability of the proposed formulation to other geological settings, lithologies, and logging configurations remains to be further validated. In addition, while the residual component is interpreted as representing structured bias, its generality across different datasets requires additional investigation.

Future work will focus on evaluating the transferability of the proposed approach across multiple basins and geological conditions, as well as extending the formulation to other heterogeneous subsurface characterization problems.

Overall, this study provides a structured approach for incorporating lithological heterogeneity into geomechanical characterization and demonstrates its effectiveness in reducing bias and improving cross-well consistency under the conditions considered.