Geological reconstruction method based on physical information neural network and anisotropic regularization

Abstract

The application discloses a geological reconstruction method based on physical information neural network and anisotropic regularization, and belongs to the technical field of geological modeling and engineering geology. The geological reconstruction method aims at a two-dimensional lithology reconstruction problem under sparse drilling conditions, maps discrete lithology into a continuous differentiable category probability field, and constructs a neural network model with spatial coordinates as input and the probability field as output. Under the maximum posterior framework, a data consistency term and an anisotropic structure regularization term are jointly optimized. In the transverse direction, a second-order smoothing is adopted to suppress high-frequency oscillation, and in the longitudinal direction, a first-order weak constraint is adopted to allow limited jumps between layers, so as to depict the characteristics of transverse continuity and longitudinal sudden change of sedimentary strata. Based on the posterior probability distribution, an entropy type uncertainty field is constructed to quantitatively evaluate the prediction confidence. The application fuses data driving and geological priori, can stably reconstruct the stratum structure under sparse drilling, maintain the layered continuity and interface rationality, and provide uncertainty information with geological significance.

Description

Technical Field

[0001] This invention relates to the fields of geological modeling and engineering geology, and in particular to a geological reconstruction method based on physical information neural networks and anisotropic regularization. Background Technology

[0002] Lithological spatial modeling is a core issue in underground engineering design, tunnel construction, and resource exploration. Due to the high heterogeneity and spatial complexity of geological bodies, practical engineering projects typically only obtain sparse observation data through a limited number of boreholes, resulting in the inversion problem of "point observations—overall unknown." How to reasonably recover the continuous stratigraphic structure under the constraint of limited boreholes has always been a key challenge in the fields of engineering geology and geoscience modeling. Traditional methods such as linear interpolation, Kriging, and spatial statistical models based on Markov random fields (MRF) can characterize certain spatial correlations, but they are prone to oversmoothing or structural fractures under sparse borehole conditions, making it difficult to simultaneously ensure data fitting and stratigraphic geometric rationality.

[0003] In recent years, machine learning and deep learning methods have been increasingly applied to geological modeling. Neural network-based classification models can improve prediction accuracy through nonlinear mapping, but most methods are purely data-driven frameworks, lacking explicit constraints on spatial structural patterns. When supervised samples are highly sparse, the model is prone to unstable extrapolations in unobserved areas, even exhibiting non-physical oscillations or isolated patches. Therefore, in geological modeling problems, simply relying on data fitting is often insufficient to guarantee structural rationality; introducing physical or geological prior constraints becomes a necessary approach. Physically-Informed Neural Networks (PINNs) embed physical laws into the loss function, achieving joint optimization of data and constraints, providing a new research approach for complex inversion problems. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a geological reconstruction method based on physical information neural networks and anisotropic regularization. This method integrates data-driven approaches with geological priors, enabling stable reconstruction of stratigraphic structures under sparse borehole conditions, maintaining layered continuity and interface rationality, and providing geologically significant uncertainty information.

[0005] To solve the above-mentioned technical problems, the technical solution proposed by this invention is as follows:

[0006] A geological reconstruction method based on physical information neural networks and anisotropic regularization includes the following steps:

[0007] Step S1: Obtain spatial discrete grids and sparse borehole observation data for a two-dimensional or three-dimensional geological domain. The sparse borehole observation data includes the borehole observation location and its corresponding lithology.

[0008] Step S2: Construct a neural network model with spatial coordinates as input and lithology category probability field as output. The output of the neural network model is mapped to a category probability vector that satisfies the probabilistic simplex constraint through the Softmax function.

[0009] Step S3, set the supervision item: data consistency item, based on the negative log-likelihood or cross-entropy loss at the borehole observation location, to fit the sparse lithology label;

[0010] Step S4, set the prior terms: anisotropic structure regularization term, apply asymmetric gradient constraints along different spatial directions in the probability field, wherein a second-order smooth constraint is used along the first direction to suppress high-frequency oscillations, and a first-order weak constraint is used along the second direction to allow finite jumps.

[0011] Step S5: Construct the joint objective function: Combine the data consistency term and the anisotropic structure regularization term to obtain the final optimization problem;

[0012] Step S6: Calculate the gradient of the neural network model output relative to spatial coordinates using automatic differentiation, then calculate the anisotropic structure regularization term, and minimize the joint loss function using stochastic gradient descent or adaptive optimization algorithms to obtain the optimal network parameters.

[0013] Step S7: Input any spatial coordinates within the geological domain into the trained neural network model to obtain a continuously differentiable class probability field, and output the lithology prediction field through the maximum a posteriori criterion.

[0014] A further improvement to the above technical solution is as follows:

[0015] Preferably, the anisotropic structure regularization term adopts a continuous functional form, which includes two terms: the first term is a second-order transverse regularization, which integrates the square of the transverse gradient and multiplies it by the transverse weight to suppress high-frequency transverse oscillations; the second term is a first-order longitudinal regularization, which integrates the absolute value of the longitudinal gradient and multiplies it by the longitudinal weight to allow finite jumps at the layer interface; wherein the transverse weight is greater than the longitudinal weight.

[0016] Preferably, the anisotropic structure regularization term is numerically approximated in discrete computation by averaging the global set of sampling points, and the gradient term is obtained by directly differentiating the probability field output by the network through automatic differentiation, without the need to construct an explicit finite difference operator.

[0017] Preferably, the neural network model is implemented using a multilayer perceptron, and the activation function of the hidden layer is the hyperbolic tangent function.

[0018] Preferably, the joint objective function is understood from the perspective of Bayesian maximum a posteriori estimation as minimizing the sum of the negative log-likelihood and the negative log-prior, wherein the prior assumption is an exponential form related to the regularization term of the anisotropic structure, and the regularization weight is equivalent to the prior strength parameter.

[0019] Preferably, the regularization weights in the joint objective function are used to control the trade-off between data fitting and structural smoothing. When the regularization weights approach 0, the model degenerates into a purely data-driven model, and when the regularization weights approach infinity, the model tends to become oversmooth.

[0020] Preferably, the step S7 of outputting the lithology prediction field using the maximum a posteriori criterion specifically involves: taking the category index corresponding to the highest probability, and then obtaining the final lithology prediction through the mapping function from the category index to the original lithology label.

[0021] Preferably, the method further includes: constructing an entropy-type uncertainty field based on the complete posterior probability distribution output by the trained neural network model, wherein the entropy-type uncertainty field is defined as the sum of negative class probabilities multiplied by the logarithms of the probabilities, and a numerical stabilization term is added to represent the prediction uncertainty of spatial location.

[0022] Preferably, the entropy-type uncertainty field has the following properties: when the probability of a certain category is close to 1, the uncertainty approaches 0; when the probabilities of multiple categories are close to a uniform distribution, the uncertainty approaches the logarithm of the number of categories; the larger the entropy value, the higher the prediction uncertainty.

[0023] The geological reconstruction method based on physical information neural networks and anisotropic regularization provided by this invention has the following advantages compared with existing technologies:

[0024] (1) The geological reconstruction method based on physical information neural network and anisotropic regularization of the present invention proposes a modeling framework of physical information neural network (PINN) based on probability field parameterization for the two-dimensional lithology reconstruction problem under sparse borehole conditions. By transforming discrete lithology mapping into a continuously differentiable class probability field, and jointly optimizing the data consistency term and the anisotropic structure regularization term under the maximum a posteriori (MAP) framework, a unified expression of data constraints and geological structure priors is achieved. The constructed transverse second-order smoothing and longitudinal first-order weak constraint mechanism effectively characterizes the structural features of sedimentary strata that are "transversely continuous and longitudinally abruptly changeable".

[0025] (2) The geological reconstruction method based on physical information neural networks and anisotropic regularization of the present invention, based on the complete posterior probability distribution of the output, further constructs an entropy-type uncertainty field, realizing a quantitative assessment of prediction confidence. High uncertainty values ​​are mainly distributed at interlayer interfaces and sparse borehole regions, possessing clear geological interpretation significance. The method of the present invention provides a unified framework for geological modeling that integrates data-driven and physical priors, offering certain advantages in terms of theoretical interpretability and engineering applicability. Attached Figure Description

[0026] Figure 1 This is a flowchart of the present invention.

[0027] Figure 2 The geological profile is shown in the experimental verification diagram, where (a) is the original strata and (b) is the known borehole.

[0028] Figure 3 The diagram shows the lithological field prediction results using the method of this invention in the experimental verification, where (a) is the MAP and (b) is the uncertainty field.

[0029] Figure 4 The results are used to verify the predictions of other models used in the experiment, where (a) is the MRF model and (b) is the IC-XGBoost model. Detailed Implementation

[0030] The following provides a detailed description of specific embodiments of the present invention. It should be understood that the specific embodiments described herein are for illustrative and explanatory purposes only and are not intended to limit the scope of the invention.

[0031] The geological reconstruction method based on physical information neural network and anisotropic regularization of the present invention represents discrete lithology as a continuously differentiable class probability field, and constructs data consistency term and anisotropic structure regularization term under the maximum a posteriori (MAP) framework to perform directional constraints on sedimentary layer structure characteristics, thereby realizing the spatial reconstruction of geological profile lithology under sparse borehole conditions.

[0032] Specifically, the following steps are included:

[0033] Step S1: Obtain spatial discrete grids and sparse borehole observation data for a two-dimensional or three-dimensional geological domain. The sparse borehole observation data includes the borehole observation locations and their corresponding lithological categories.

[0034] Specifically, it includes the following:

[0035] (1) Spatial discrete grid.

[0036] Considering the spatial domain of two-dimensional geological profiles ,in This represents two-dimensional Euclidean space. Spatial position is determined by coordinate vectors. This indicates that, to eliminate the impact of scale differences on network training, spatial coordinates are normalized to a unit interval. ,in Indicates the horizontal direction of the cross-section. Indicates the vertical (depth) direction.

[0037] In discrete implementation, the research domain is divided into The regular grid, corresponding to the coordinate set , , It is a discrete grid. Wherein, This represents the number of vertical discrete nodes (number of grid rows). The number of discrete nodes (grid columns) in the horizontal direction, and the total number of grid points is: . Indicates the first Normalized horizontal coordinates of each grid point; Indicates the first Normalized vertical coordinates of grid points.

[0038] (2) Sparse borehole observation data.

[0039] Lithology is a discrete categorical variable, and its values ​​belong to a finite set. , This represents the c-th lithological category (c=1, 2, 3, ...). ), This represents the total number of lithological categories. Therefore, the lithological field can be considered as defined in a spatial domain. Discrete mappings on Because lithological labels are only observable at borehole locations, complete discrete mapping is not possible. Since it is unknown, it needs to be rebuilt through learning.

[0040] Borehole observations are only available at a limited number of locations, forming a set of observation points. , , .in This is the location for borehole observation. Lithological tags; The total number of observed samples, The sample index is (n = 1, 2, …, N), and the number of samples is much smaller than the number of grid points in the global domain. Unobserved locations are not included in the calculation of supervised loss, but will be used as objects for applying physical constraints later to enhance the global generalization ability.

[0041] Step S2: Construct a neural network model with spatial coordinates as input and lithology category probability field as output. The output of the neural network model is mapped to a category probability vector that satisfies the probabilistic simplex constraint through the Softmax function.

[0042] Specifically, the following steps are included:

[0043] S2-1, Lithology category probability field.

[0044] Direct learning of discrete mappings The lack of differentiable structure hinders the introduction of physical constraints. Therefore, a class probability field representation is introduced:

[0045] (1)

[0046] in, Indicates spatial location The category probability vector at that location, i.e., the probability field; For category exist The probability at a given location satisfies and ; The total number of lithological categories, i.e., the category set; where , For probability operators, This represents the i-th lithological category. Spatial location. Lithology predictions can be obtained using the maximum a posteriori (MAP) criterion.

[0047] S2-2, Construct a neural network model.

[0048] In the continuous spatial domain By constructing a differentiable class probability field, this invention uses a neural network based on coordinate input to parameterize the probability distribution, transforming the discrete lithology inversion problem into a parameter optimization problem in a continuous function space.

[0049] S2-2-1, Coordinate-driven function mapping.

[0050] Define parameters Neural network mapping:

[0051] (2)

[0052] in, Indicates Neural network functions with parameters It is a C-dimensional real vector space, consistent with the lithology category. The unnormalized score (logits) output by the network and , This represents the unnormalized score (logit) of category c at position x. Represents the output vector of the neural network at point x; parameters Represents the complete set of learnable parameters of the network, where and These are the weight matrix and the bias vector, respectively. Indicates the total number of network layers. This represents the hidden layer. This mapping essentially defines a family of continuous functions to describe the relative confidence of each category in the space.

[0053] S2-2-2, Probabilistic Simplex Mapping (Softmax Probabilistic Construction).

[0054] To ensure the output satisfies the probability constraints, logits are mapped to a probabilistic simplex using a softmax transformation. For the th... Classes include:

[0055] (3)

[0056] in, Let Z represent the conditional probability under parameter θ, and let Z denote the random variable (lithology category). It is the i-th type of lithology. For category indexing, For the logit of category i, The logit for category j.

[0057] Therefore, the parameterized class probability field can be defined as:

[0058] (4)

[0059] in, Indicated by neural network parameters The parameterized class probability field, Output probabilities for the model. for The probabilistic simplex is expressed as:

[0060] (5)

[0061] in, It is a probability vector. Let i be the i-th component of the probability vector.

[0062] This parameterization method has four properties: (1) The output naturally satisfies non-negativity and normalization; (2) The probability field is continuously differentiable with respect to spatial coordinates; (3) It is suitable for multi-class discrete field modeling; (4) It is convenient to introduce gradient-type physical regularization terms later.

[0063] Therefore, the parameterized class probability field can be expressed as:

[0064] (6)

[0065] in, This represents the Softmax operator.

[0066] S2-2-3, network structure form.

[0067] function It is implemented using a multilayer perceptron (MLP). For the th Hidden layers:

[0068] (7)

[0069] in, For the first Hidden representation (feature vector) of a layer For the first The weight matrix of the layer, For the first Hidden representation (feature vector) of a layer For the first Layer bias vector, ; This is the activation function.

[0070] The output layer is a linear mapping:

[0071] (8)

[0072] in, This indicates the unnormalized score. The weight matrix of the output layer. Output for the last hidden layer. This is the bias vector for the output layer. This represents the total number of network layers (including the output layer).

[0073] This invention selects As the activation function of the hidden layer. Because It provides a smooth, continuously differentiable nonlinear mapping, and in the PINN framework, it is beneficial for the stable computation of spatial gradients and higher-order derivatives, and also helps to reduce gradient oscillations and numerical instability problems.

[0074] Step S3, set the supervision item: data consistency item, based on the negative log-likelihood or cross-entropy loss at the borehole observation location, to fit the sparse lithology label.

[0075] S3-1, Statistical Modeling Perspective.

[0076] In the set of observation points Above, assuming lithological labels At a given spatial borehole observation location and model parameters Under the condition that it follows a multinomial distribution, its class probability is determined by the network parameterized class probability field. Given, that is:

[0077] (9)

[0078] in, Represents a probability operator. In model parameters Next position Category The probability of.

[0079] Under the conditional independence assumptions of spatial location and model parameters, the likelihood function of the observed data is:

[0080] (10)

[0081] in, For the observed data about the parameters The likelihood function, For the model at position For the true category The predicted probability.

[0082] To perform parameter estimation, the maximum likelihood (MLE) principle is adopted, which is equivalent to minimizing the negative log-likelihood (NLL):

[0083] (11)

[0084] in, This represents the data consistency loss, which is the negative log-likelihood of the observed data under the model parameters. It is used to measure the difference between the predicted distribution and the true label. This form is equivalent to the multi-class cross-entropy loss.

[0085] Within the inversion framework, This constitutes the data fidelity term in the objective functional. This data fidelity term increases the posterior probability of the true class. This causes the parameterized lithology probability field to approximate the true label distribution at the observation points. Specifically, it applies stronger constraints at the borehole locations to maximize the predicted probability at the true category. From a Bayesian perspective, this data consistency term corresponds to the observed data in the model parameters. negative log-likelihood That is, the maximum likelihood estimation (MLE) objective.

[0086] S3-2, Limitations of sparse supervision.

[0087] In practical geological problems This means that borehole observations are highly sparse relative to the global grid. In this case, relying solely on... Estimating model parameters solely by fitting borehole observation data often leads to the following problems: (1) instability in extrapolating unobserved areas; (2) non-physical high-frequency oscillations in the probability field; (3) discontinuous lithological spatial structure; and (4) the appearance of isolated patches or noise patterns. In other words, It only guarantees the consistency of point states, but cannot constrain the rationality of spatial structure.

[0088] Therefore, in the context of sparse data, this invention introduces structural priors or physical constraints to achieve joint optimization of data consistency and structural rationality.

[0089] S3-3, Function Space.

[0090] At the function space level, the above supervision terms can be viewed as discrete sampling constraints on the probability field:

[0091] (12)

[0092] in, This is the cross-entropy loss function. Therefore, the data items are not directly restricted. Throughout The behavior is constrained only on a finite number of sample points. This property is the fundamental reason why subsequent physical regularization terms are needed.

[0093] Step S4, set physical / geological priors: anisotropic structure regularization terms, and apply asymmetric gradient constraints in different spatial directions on the probability field, wherein a second-order smooth constraint is used along the first direction (lateral) to suppress high-frequency oscillations, and a first-order weak constraint is used along the second direction (longitudinal) to allow finite jumps.

[0094] S4-1, Structural a priori motivation.

[0095] Sedimentary strata typically exhibit an approximately layered structure in space, along the strike (horizontal) The direction has strong continuity; in depth (vertical) Relative abrupt changes can occur in the direction (interlayer interface).

[0096] The spatial structure of lithology satisfies the anisotropic smoothness assumption:

[0097] (13)

[0098] Since the lithological field has been converted into a probability field This structural knowledge should be reflected in the directional differences of the probability field gradient.

[0099] Unlike traditional PINN, which constructs strong physical constraints by minimizing the residuals of the control partial differential equations (PDEs), the physical information introduced in this invention originates from the statistical structural regularities of sedimentary strata and can be regarded as a soft constraint for geological consistency. Its goal is to avoid generating non-physical lateral oscillations that violate layered continuity while ensuring borehole data consistency.

[0100] S4-2, continuous functional form.

[0101] Due to probability field Since spatial coordinates are continuously differentiable, anisotropic gradient functionals can be constructed:

[0102] (14)

[0103] in, It is an anisotropic structure regular functional used to impose directional constraints on the spatial gradient, thereby introducing the prior knowledge of the stratigraphic structure during the model optimization process; These are the lateral gradient weights, controlling the intensity of lateral smoothing. These are the vertical regularization weights, controlling the strength of the vertical constraints; they typically satisfy... This indicates stronger lateral constraints.

[0104] The functional has the following properties: (1) The first term is a transverse second-order (Tikhonov-type) regularization, which suppresses transverse high-frequency oscillations; (2) The second term is a longitudinal first-order (TV-type) weak regularization, which allows for finite jumps at the layer interface; (3) Both terms together characterize the directional differences of the deposited layered structure. Therefore, the regularization can be regarded as an anisotropic probability field smoothing functional.

[0105] S4-3, Discrete Approximation and Numerical Implementation.

[0106] Because the actual calculation is performed on a finite set of coordinates The continuous functional can be discretized using Monte Carlo averaging:

[0107] (15)

[0108] in, For local energy density, Indicates spatial sampling distribution The mathematical expectation is the average of the sampled points over the spatial domain Ω.

[0109] The expression for local energy density is:

[0110] (16)

[0111] Therefore, the numerical implementation can be obtained as follows:

[0112] (17)

[0113] The gradient term is directly applied to the network output through automatic differentiation. Differentiation is obtained without constructing an explicit finite difference operator, thus avoiding truncation errors and numerical noise introduced by grid discretization. Furthermore, although the category dimension is aggregated at the implementation level (i.e., for...),... (Summation), the essence of this regularization term is still to impose asymmetric constraints on the rate of change of the probability field in different spatial directions. Second-order smoothing is used in the transverse direction to suppress high-frequency oscillations, while first-order weak constraints are used in the longitudinal direction to allow the existence of interlayer interfaces, thereby characterizing the anisotropic features of sedimentary stratigraphic structures.

[0114] S4-4, The rationale for probability space regularization.

[0115] This invention applies structure regularization to the probability space rather than the logits space. If directly in the logits space... Constructing gradient constraints on the upper part may lead to two problems due to the nonlinear characteristics of the Softmax transformation: (1) The class competition relationship is distorted: the smoothing of logits is not equivalent to the smoothing of the posterior probability, which may destroy the relative proportion structure between classes; (2) Gradient amplification effect: the exponential mapping of Softmax may nonlinearly amplify the small perturbations of logits, thereby causing unstable gradient propagation in the probability space.

[0116] In contrast, this invention is in the probability field Applying regularity to the above has a clearer statistical and physical meaning, which is reflected in: (1) the probability is naturally bounded ( (1) Numerical scale stability; (2) The gradient directly corresponds to the spatial rate of change of the posterior distribution; (3) The object of the constraint is the final lithological probability field, and the physical interpretation is more intuitive.

[0117] Therefore, the anisotropic smoothing regularization is directly applied to It is mathematically more stable, statistically more consistent, and physically more in line with the modeling objectives of lithological probability fields.

[0118] Step S5: Construct the joint objective function.

[0119] Combining the data consistency term and the anisotropic structure regularization term, the final optimization problem of this invention can be uniformly expressed as:

[0120] (18)

[0121] in, For the optimization objective of the model parameter θ, The negative log-likelihood corresponds to the borehole observation constraint. As the regularization weight, and Used to control the trade-off between data fitting and structural smoothing.

[0122] The objective function can be understood from the perspective of Bayesian maximum a posteriori estimation (MAP):

[0123] (19)

[0124] in, Indicates the optimal model parameters. The logarithmic form of the prior distribution of parameters reflects the model's prior assumptions about structural smoothness and spatial continuity.

[0125] The prior assumption is:

[0126] (20)

[0127] Therefore, the regularization weight λ is equivalent to the prior strength parameter. When When the size is large, structural priors dominate; when When the size is small, the data items dominate. The entire model framework is essentially a MAP probability field inversion model under sparse supervision.

[0128] The regularization weight λ controls the trade-off between fitting error and structural complexity, which can be understood as a bias-variance trade-off mechanism: λ→0: degenerates into a purely data-driven model; λ→∞: tends to oversmooth.

[0129] Step S6: Calculate the gradient of the neural network model output relative to spatial coordinates using automatic differentiation, then calculate the anisotropic structure regularization term, and minimize the joint loss function using stochastic gradient descent or adaptive optimization algorithms to obtain the optimal network parameters.

[0130] The continuous variational optimization problem is implemented within a discrete sampling framework, specifically as follows:

[0131] S6-1, global sampling and automatic differentiation are implemented.

[0132] because For a function that is continuously differentiable with respect to spatial coordinates, the spatial gradient term involved in anisotropic canonical functionals is... , It can be directly calculated through automatic differentiation. Unlike the traditional finite difference method, this invention does not construct an explicit difference operator, but automatically differentiates the network output. This implementation has the following advantages: (1) avoiding the truncation error caused by grid difference; (2) eliminating the noise amplification effect introduced by numerical difference; (3) maintaining the continuous differentiable structure of the probability field; (4) realizing the natural coupling of regular gradient and parameter update.

[0133] In discrete computation, the continuous integral form (Equation (14)) is numerically approximated by averaging over the entire set of sampled points (step S4-3). Therefore, the continuous variational problem is transformed into an energy minimization problem based on discrete sampling.

[0134] S6-2, Optimize strategies and training process.

[0135] The joint objective function of the final optimization problem includes both strong local constraints and soft global constraints, with the strong local constraints being the data consistency term. Global soft constraints are anisotropic structural regularization. .

[0136] Parameter optimization is performed using the first-order stochastic gradient method (Adam optimizer). The learning rate controls the magnitude of parameter updates in each iteration.

[0137] Each training iteration includes the following steps:

[0138] (1) Forward propagation of borehole sample points: Calculate the data consistency loss according to equation (11);

[0139] (2) Forward propagation of global sampling points: Calculate the probability field according to equation (6). ;

[0140] (3) Automatic differentiation to calculate spatial gradient: used to construct the anisotropic energy term in equation (14);

[0141] (4) Construct the total loss function: Combine the data terms and regularization terms according to equation (18);

[0142] (5) Backpropagation and parameter update.

[0143] This process achieves joint optimization of point state supervision (drilling constraints) and spatial structure constraints (sedimentary layered priors).

[0144] It is important to emphasize that this invention applies structural regularization to the probability space rather than the logits space. Therefore, the regularization term directly acts on the probability field defined by Equation (6); the probability values ​​are naturally bounded (0~1), and the numerical scale is stable; the gradient corresponds to the spatial rate of change of the posterior distribution, and the physical interpretation is clear; the problem of Softmax nonlinearly amplifying logits perturbations is avoided. Within the entire PINN model framework, the optimization process of Equation (18) can be regarded as a direct control of the family of probability field functions, thereby achieving the unity of data consistency (Equation (11)) and structural rationality (Equation (14)).

[0145] Step S7: Input any spatial coordinates within the geological domain into the trained neural network model to obtain a continuously differentiable class probability field, and output the lithology prediction field through the maximum a posteriori criterion.

[0146] This invention formulates the lithology inversion problem as a minimization problem of a joint objective function, provides a theoretical explanation from the perspective of MAP, and further illustrates the numerical solution process of this variational problem.

[0147] The parameter optimization converges and the optimal solution is obtained. Then, the model defines a space domain. A continuously differentiable probability field:

[0148] (twenty one)

[0149] in, Indicates the optimal parameters The category probability vector is the posterior probability distribution of the lithology category at spatial location x. This represents the probability of category c at position 𝐱. for A 3D probability simplex. This probability field is the optimal posterior distribution approximation in the sense of MAP.

[0150] Based on this probability field, lithology prediction results and their uncertainty measures can be further obtained.

[0151] S7-1 Maximum a posteriori lithology field (MAP lithology field).

[0152] Since the model outputs the complete class probability distribution, rather than a single label, obtaining the optimal parameters... Then, for any spatial location Lithological prediction can be determined using the maximum a posteriori (MAP) criterion:

[0153] (twenty two)

[0154] in, This is the lithology prediction result at location x. The category index corresponding to the highest probability; A mapping function for category indexes to original lithology labels.

[0155] In discrete grid By reasoning from the above, a two-dimensional lithology prediction field can be obtained. . Indicates that it is defined in A category field on a grid, where the value at each position belongs to a category set. .

[0156] The predicted field is the optimal reconstruction result of the joint optimization problem in the sense of MAP.

[0157] From the perspective of function space, this result is a pointwise decision mapping of a continuous probability field, thus maintaining a balance between data consistency and structural rationality.

[0158] S7-2, Uncertainty Quantification Based on Entropy.

[0159] Unlike traditional deterministic interpolation methods, the model in this invention outputs a complete posterior class probability distribution. Therefore, while obtaining the MAP prediction, the prediction confidence can also be quantified.

[0160] This invention uses Shannon information entropy as an indicator of classification uncertainty:

[0161] (twenty three)

[0162] in, This is the prediction uncertainty function defined based on Shannon's information entropy. It is a numerically stable term; It indicates spatial location The prediction uncertainty has the following properties: (1) When the probability of a certain category is close to 1, This indicates that the prediction is highly certain; (2) when the probabilities of multiple categories are close to a uniform distribution: This indicates that the model has difficulty distinguishing lithological categories. Therefore, the higher the entropy value, the higher the prediction uncertainty. The uncertainty field can be obtained on a discrete grid. , for The space of real-valued vectors. This uncertainty field originates from the posterior probability distribution of the categories themselves and belongs to predictive uncertainty, rather than the uncertainty of the parametric posterior distribution.

[0163] From a geological perspective, uncertainties typically possess clear spatial structural significance: near inter-stratum interfaces; lithological transition zones; regions with sparse borehole information; and areas of significant multi-class probability competition. Because this invention applies anisotropic structural regularization to the probability space, high uncertainty regions often correspond to structural transition zones or regions with weaker data constraints. Therefore, the uncertainty results not only provide model confidence information but can also be used to: assist in identifying potential stratigraphic interfaces; guide subsequent borehole placement; evaluate the reliability of model extrapolation; and analyze the balance between structural priors and data constraints.

[0164] From the perspective of MAP interpretation (Equation (19)), the entropy field reflects the degree of uncertainty of the optimal posterior distribution in the prediction space, and is a natural output of the probability field modeling framework.

[0165] Experimental verification:

[0166] To verify the applicability and effectiveness of the probabilistic field geological reconstruction model based on Physics-Informed Neural Network (PINN) proposed in this invention in practical engineering, a geological profile example of a tunnel project in Australia reported in Shiand Wang (2021) was selected as the research object. In this case, boreholes were arranged along the tunnel axis for geological exploration. The strata were mainly controlled by sedimentation and erosion, exhibiting obvious layered structure characteristics, but local strata showed wedge pinch-outs and interlayer cross-linking.

[0167] The overall dimensions of the geological profile are 600m in horizontal length and 40m in vertical depth. For example... Figure 2 As shown in (a), the site mainly comprises four stratigraphic types (from top to bottom): the Bassendean Sand (BS), the Guildford Formation (GF), the Gnangara Sand (GS), and the Ascot Formation (AF). BS and GS are primarily sand layers, while GF and AF are relatively dense strata. GS is absent in some boreholes, creating stratigraphic discontinuities. The complex morphology of the GS / GF and GS / AF interfaces presents a challenge for profile interpretation. This phenomenon of a layer appearing / disappearing in adjacent boreholes is a typical problem that traditional straight-line connection methods struggle to explain reasonably.

[0168] To facilitate numerical computation and model training, the original geological profile was discretized using a regular grid: a horizontal resolution of 3m and a vertical resolution of 0.4m, resulting in a grid size of 100×50. In the case study verification phase, only four equally spaced boreholes were selected, such as... Figure 2As shown in (b), the borehole spacing is 100m, and the borehole data accounts for about 4% of all grid points, which forms a typical problem of sparse observation and global unknown reconstruction.

[0169] Using the PINN model framework constructed in this invention, the Australian tunnel profile case was trained, and the lithological prediction results and uncertainty distribution were obtained by reasoning on the global grid.

[0170] Regarding MAP and uncertainty quantification results: Figure 3 (a) The model-predicted MAP lithological field results are presented. It can be observed that, using only four boreholes (approximately 4% of the observations), the model successfully recovered the overall layered structure characteristics, and the lithological distribution is highly consistent with the actual profile morphology. Quantitative evaluation results show that the overall classification accuracy reaches Acc=95.16%, indicating that the proposed probabilistic field modeling and anisotropic regularization mechanism can achieve stable extrapolation under sparse supervision.

[0171] From the perspective of spatial structure, the prediction results exhibit the following characteristics:

[0172] (1) Good lateral continuity

[0173] The strata showed obvious continuity along the horizontal direction, and no high-frequency oscillations or isolated patches were observed, indicating that the second-order smoothing term in the horizontal direction (the first term of Equation (14)) effectively suppressed non-physical oscillations.

[0174] (2) The interface form is reasonably preserved.

[0175] At the junction of GS and AF, the model successfully depicts the wedge-shaped interface structure sloping downwards to the right, without being overly smoothed. This indicates that the longitudinal first-order weak regularization term allows for limited jumps in the layer interface, preserving the characteristics of the depositional structure.

[0176] (3) Complex transition zones can be reconstructed

[0177] The GS layer shows a tendency to thin or even disappear in local areas, and the prediction results can capture this pinch-out phenomenon, demonstrating the model's sensitivity to interlayer changes.

[0178] It is worth noting that the actual boundary shown by the dashed line largely coincides with the predicted interface, with only a slight offset in the rapidly changing interlayer regions. This error is mainly concentrated in sections with large borehole spacing, consistent with the extrapolation error distribution under sparse monitoring conditions.

[0179] Figure 3(b) presents the corresponding uncertainty field. It can be observed that the uncertainty distribution exhibits distinct structural characteristics: high uncertainty regions are mainly concentrated near interlayer interfaces; uncertainty increases significantly in the GS pinch-out region; and local probability competition exists in areas with large borehole spacing. When the probability of a certain category approaches 1, the entropy value approaches 0; when the probabilities of multiple categories are close to a uniform distribution, the entropy value increases. Therefore, the entropy field reflects the dispersion of the predicted posterior distribution, indicating that the uncertainty quantification results have clear geological significance.

[0180] Regarding the verification of the structural prior effects: combined with the anisotropic regularity, it can be seen that no oscillating fringes appear in the transverse direction; abrupt changes in the interface are allowed in the longitudinal direction; the probability field is smooth but not excessively blurred.

[0181] This demonstrates that the joint optimization mechanism of data consistency and prior geological structure effectively achieves a balance between data-driven approaches and physical constraints.

[0182] To further verify the effectiveness of the PINN probability field model proposed in this invention, it is compared and analyzed with two representative data-driven methods: one is a spatial statistical model based on Markov random fields (MRF), and the other is an IC-XGBoost model based on feature learning. The three methods were trained and predicted under the same drilling conditions, and the results are as follows: Figure 4 As shown.

[0183] (1) Overall accuracy comparison

[0184] In terms of overall classification accuracy, the model of this invention achieves 95.16%, the MRF model 94.8%, and the IC-XGBoost model 91.2%. It can be seen that, under the same sparse supervision conditions, the method of this invention is slightly better than the MRF model in overall prediction accuracy and significantly better than the IC-XGBoost model.

[0185] (2) Contrast of layered continuity

[0186] Regarding lateral continuity, the MRF model can maintain the layered structure to some extent through neighborhood correlation constraints, but block splicing still occurs in sections with large borehole spacing, and the interface shows slight step-like fluctuations locally. The IC-XGBoost model, due to the lack of explicit spatial structure constraints, produces obvious class jumps between boreholes, has weak layer continuity, and exhibits non-physical fractures in some areas.

[0187] In contrast, the model of this invention exhibits a smooth, continuous layered structure across the entire cross-section. The absence of high-frequency oscillations or isolated patches along the horizontal direction indicates that the transverse second-order smoothing mechanism in the anisotropic regularization term effectively suppresses non-physical oscillations. This continuity is highly consistent with the actual characteristics of sedimentary strata.

[0188] (3) Ability to depict interlayer interfaces

[0189] In the GS–AF interface region, the actual strata exhibit obvious tilting and bending characteristics. The MRF model can roughly recover the interface trend, but local misalignment or interface blurring occurs. The IC-XGBoost model shows significant interface shift around 200m, with layer thickness distribution differing considerably from the actual structure.

[0190] The model of this invention can accurately recover interface tilt angle changes while maintaining overall interface continuity. The interface is neither overly smoothed nor exhibits abrupt, discontinuous transitions. This result indicates that the longitudinal first-order weak regularization term allows for finite abrupt changes in the interlayer interface, while the transverse smoothing term ensures overall interface stability, thus achieving a continuous yet flexible structural recovery effect.

[0191] (4) Capability to handle pinch-out and missing layers

[0192] On the right side of the profile, the GS layer gradually thins and eventually pinches out. This region is a key area where traditional methods are prone to misjudgment. The IC-XGBoost model tends to forcibly extend the GS layer to the right end, resulting in excessive extrapolation of the layer volume. The MRF model exhibits some smoothing residue in the pinch-out region and fails to fully capture the layer volume disappearance trend.

[0193] The model of this invention exhibits a gradual thinning and eventual disappearance trend in this region, which is basically consistent with the real structure, and no isolated patches or anomalous segments are generated. This shows that structural regularization in the probability space can naturally achieve layer convergence through a class competition mechanism without the need for artificial rules.

[0194] (5) Analysis of differences in modeling mechanisms

[0195] The differences among the three methods essentially stem from their different modeling paradigms.

[0196] IC-XGBoost is a purely data-driven model whose predictions rely entirely on local feature mappings and lack global constraints on spatial continuity. Therefore, it is prone to structural distortion under sparse borehole conditions.

[0197] The MRF model characterizes spatial correlation through neighborhood statistical relationships, but its smoothing mechanism is usually isotropic, making it difficult to distinguish between horizontal continuity and vertical abrupt changes, thus limiting its ability to characterize interfaces.

[0198] The model in this invention is based on continuous probability field parameterization and jointly optimizes the data consistency term and the anisotropic structure regularization term within the MAP framework. Second-order smoothing is used laterally to suppress oscillations, while first-order weak constraints are employed vertically to allow interface jumps, which mechanistically better reflects the physical characteristics of sedimentary strata. This directional structural constraint allows the model to maintain the rationality of stratigraphic geometry even under sparse supervision.

[0199] The above embodiments are merely preferred examples of the present invention and are not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Therefore, any simple modifications, equivalent changes, and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention should fall within the protection scope of the present invention.

Claims

1. A geological reconstruction method based on physical information neural networks and anisotropic regularization, characterized in that, Includes the following steps: Step S1: Obtain spatial discrete grids and sparse borehole observation data for a two-dimensional or three-dimensional geological domain. The sparse borehole observation data includes the borehole observation location and its corresponding lithology. Step S2: Construct a neural network model with spatial coordinates as input and lithology category probability field as output. The output of the neural network model is mapped to a category probability vector that satisfies the probabilistic simplex constraint through the Softmax function. Step S3, set the supervision item: data consistency item, based on the negative log-likelihood or cross-entropy loss at the borehole observation location, to fit the sparse lithology label; Step S4, set the prior terms: anisotropic structure regularization term, apply asymmetric gradient constraints along different spatial directions in the probability field, wherein a second-order smooth constraint is used along the first direction to suppress high-frequency oscillations, and a first-order weak constraint is used along the second direction to allow finite jumps. Step S5: Construct the joint objective function: Combine the data consistency term and the anisotropic structure regularization term to obtain the final optimization problem; Step S6: Calculate the gradient of the neural network model output relative to spatial coordinates using automatic differentiation technology, then calculate the anisotropic structure regularization term, and minimize the joint loss function using stochastic gradient descent or adaptive optimization algorithm to obtain the optimal network parameters. Step S7: Input any spatial coordinates within the geological domain into the trained neural network model to obtain a continuously differentiable class probability field, and output the lithology prediction field through the maximum a posteriori criterion.

2. The geological reconstruction method based on physical information neural network and anisotropic regularization according to claim 1, characterized in that, The anisotropic structure regularization term adopts a continuous functional form, which contains two terms: the first term is a second-order transverse regularization, which integrates the square of the transverse gradient and multiplies it by the transverse weight to suppress high-frequency transverse oscillations; the second term is a first-order longitudinal regularization, which integrates the absolute value of the longitudinal gradient and multiplies it by the longitudinal weight to allow finite jumps at the layer interface; wherein the transverse weight is greater than the longitudinal weight.

3. The geological reconstruction method based on physical information neural network and anisotropic regularization according to claim 2, characterized in that, The anisotropic structure regularization term is numerically approximated in discrete computation by averaging the global set of sampling points, and the gradient term is obtained by automatically differentiating the probability field output by the network without constructing an explicit finite difference operator.

4. The geological reconstruction method based on physical information neural network and anisotropic regularization according to claim 3, characterized in that, The neural network model is implemented using a multilayer perceptron, and the activation function of the hidden layer is the hyperbolic tangent function.

5. The geological reconstruction method based on physical information neural network and anisotropic regularization according to claim 4, characterized in that, The joint objective function, from the perspective of Bayesian maximum a posteriori estimation, is the minimization of the sum of the negative log-likelihood and the negative log-prior, where the prior assumption is an exponential form related to the regularization term of the anisotropic structure, and the regularization weight is equivalent to the prior strength parameter.

6. The geological reconstruction method based on physical information neural network and anisotropic regularization according to claim 5, characterized in that, The regularization weights in the joint objective function are used to control the trade-off between data fitting and structural smoothing. When the regularization weights approach 0, the model degenerates into a purely data-driven model, and when the regularization weights approach infinity, the model tends to become oversmooth.

7. The geological reconstruction method based on physical information neural network and anisotropic regularization according to claim 1, characterized in that, The step S7, which outputs the lithology prediction field using the maximum a posteriori criterion, specifically involves taking the category index corresponding to the highest probability and then obtaining the final lithology prediction through the mapping function from the category index to the original lithology label.

8. The geological reconstruction method based on physical information neural network and anisotropic regularization according to claim 1, characterized in that, The method further includes: constructing an entropy-type uncertainty field based on the complete posterior probability distribution output by the trained neural network model. The entropy-type uncertainty field is defined as the sum of negative class probabilities multiplied by the logarithms of the probabilities, and a numerical stabilization term is added to represent the prediction uncertainty of spatial location.

9. The geological reconstruction method based on physical information neural network and anisotropic regularization according to claim 8, characterized in that, The entropy-type uncertainty field has the following properties: when the probability of a certain category is close to 1, the uncertainty approaches 0; when the probabilities of multiple categories are close to a uniform distribution, the uncertainty approaches the logarithm of the number of categories; the larger the entropy value, the higher the prediction uncertainty.

Images