A geological implicit modeling and multi-source data dynamic interpolation method and system

Abstract

The application provides a geological implicit modeling and multi-source data dynamic interpolation method and system, comprising: constructing an engineering geology knowledge graph, analyzing semantic rules and establishing a "semantic mathematics" mapping, converting sequence, morphology and structural rules into boundary conditions, anisotropy parameters and penalty functions of an implicit potential field function; fusing drilling hard data and geophysical soft data evaluated by confidence based on Bayesian inference; constructing an energy functional containing data fitting, smoothing and knowledge violation penalty and iteratively solving to generate a three-dimensional geological model; and realizing local dynamic updating by using a hierarchical spatial index. The method is suitable for water conservancy and hydropower engineering, traffic tunnel and deep mine and other engineering scenes requiring high-precision three-dimensional geological modeling and dynamic updating during design and construction.

Description

Technical Field

[0001] This application relates to the interdisciplinary fields of geological engineering digitization, three-dimensional geological modeling and artificial intelligence, specifically to a method and system for implicit geological modeling and dynamic interpolation of multi-source data. Background Technology

[0002] In the construction of "digital twins" in fields such as water conservancy and hydropower projects, transportation tunnels, and deep-earth mines, high-precision 3D geological models are the core foundational data. Currently, the mainstream modeling technology in the industry is transitioning from "explicit modeling" (manual profile delineation) to "implicit modeling" (spatial interpolation based on mathematical functions). Although implicit modeling techniques (such as those based on radial basis functions (RBF) and Kriging interpolation) have improved automation, they still have the following core pain points: 1. Interpolation based solely on geometric distance leads to a lack of geological logic: Existing algorithms interpolate purely based on geometric distance, ignoring geological constraints. In sparsely drilled areas, algorithms, in pursuit of mathematical "smoothness," often generate morphologies that violate geological common sense. For example, they interpolate lenses that should extend in layers under the control of sedimentary rhythms as spherical shapes, or they forcibly interpolate continuously on both sides of a fault, ignoring stratigraphic movement. Such models, which disregard geological constraints, are difficult to use for engineering decision-making.

[0003] 2. Difficulty in fusing multi-source data: Engineering sites contain both "hard data" (high-precision borehole data) and "soft data" (low-precision, multi-solution geophysical inversion models). Existing methods struggle to scientifically fuse these two types of data within a single mathematical equation. Geophysical data is often left idle or used only as background reference, failing to effectively constrain the range of inferences between boreholes.

[0004] 3. Lagging Dynamic Updates: Implicit modeling typically involves solving large systems of linear equations (matrix inversion), which has high computational complexity. Once new data is revealed during construction, it is often necessary to recalculate the entire site's tens of millions of grid cells, which is extremely time-consuming and cannot meet the requirement of "real-time dynamic updates as construction progresses" during the construction period.

[0005] Therefore, there is an urgent need for an intelligent modeling method that can transform the prior knowledge (semantics) of geological experts into mathematical constraints, scientifically screen and integrate multi-source data, and support rapid local updates. Summary of the Invention

[0006] In view of this, this application provides a geological implicit modeling and multi-source data dynamic interpolation method and system, which aims to solve the problems in the existing technology that are caused by the lack of geological logic due to "data-only" approach, the difficulty in scientifically integrating multi-source exploration data, and the inability to respond to real-time corrections during the construction period due to the lag in dynamic model updates.

[0007] This application proposes a method for implicit geological modeling and dynamic interpolation of multi-source data, characterized by the following steps: Acquire multi-source geological exploration data of the target area, and construct or call an engineering geological knowledge graph, which covers geological entities, geological attributes, and semantic relationships between entities; The geological semantic rules in the engineering geology knowledge graph are analyzed, and the mapping relationship between semantic rules and mathematical constraints is established. The sequence relationship rules, geometric morphology rules and structural topology rules are transformed into the boundary conditions, anisotropic parameters and penalty functions of the implicit potential field function, respectively. The multi-source geological exploration data is fused based on Bayesian inference theory. The borehole exploration data is defined as deterministic hard constraints, and the geophysical inversion data is transformed into probabilistic soft constraints after confidence evaluation using the applicability rules in the engineering geology knowledge graph. An implicit potential field energy functional is constructed, which includes a data fitting term, a smoothing term, and a knowledge violation penalty term. The knowledge violation penalty term is used to measure the degree to which the implicit potential field function violates the geological semantic rules. The energy functional is solved by an iterative optimization algorithm to obtain the implicit function expression of the geological interface. Based on the solved implicit function, isosurfaces are extracted to generate a three-dimensional geological model, and a hierarchical spatial index structure is constructed. The three-dimensional geological model is then dynamically updated locally based on newly added exploration data.

[0008] Furthermore, the specific method for transforming the hierarchical relationship rule into boundary conditions is as follows: If the engineering geological knowledge graph defines stratum A as overlying stratum B, then an inequality constraint is constructed in the normal direction of the stratum contact surface, which forces the implicit potential field value of stratum A at any spatial point to be greater than the potential field value of stratum B, in order to prevent non-geologically induced stratum inversion phenomena from occurring during the modeling process.

[0009] Furthermore, the specific method for converting geometric morphological rules into anisotropic parameters is as follows: Based on the geological body genetic types and occurrence attributes recorded in the engineering geological knowledge graph, including strike, dip and dip angle, an anisotropic covariance matrix or weight matrix is ​​constructed. The matrix is ​​used to control the range ratio and rotation angle of the interpolation kernel function in different spatial directions so that the generated geological body morphology conforms to its sedimentary or tectonic background. Specifically, for sedimentary geological bodies extending in a layered manner, the matrix makes the range of the interpolation kernel function along the strike-dip direction of the bedding plane significantly greater than the range along the normal direction; for intrusive geological bodies occurring in vein-like patterns, the matrix makes the interpolation kernel function flattened into an ellipsoid, consistent with the occurrence of the vein.

[0010] Furthermore, the specific method for transforming the construction topology rule into a penalty function is as follows: For fault cutting relationships, a step function or drift term related to the fault distance is introduced into the implicit equation on both sides of the fault to express the displacement discontinuity of the strata on both sides of the fault. Based on the motion properties of the fault, a trend rotation matrix is ​​added to the drift term to simulate the attitude rotation and traction deformation of the hanging wall relative to the footwall, thereby achieving a complete mathematical expression of the fault cutting relationship and associated structural deformation.

[0011] Furthermore, the specific steps for assessing the confidence level of the geophysical inversion data and converting it into soft constraints are as follows: Query the geophysical method applicability rules and the mapping rules between physical properties and lithology in the engineering geology knowledge graph, and score the confidence level of each anomaly area in the geophysical inversion data. Remove noise anomalies with a confidence level below a preset threshold; The values ​​of the retained high-confidence anomaly areas are normalized to the range of zero to one to characterize the probability of the corresponding geological attributes existing in that spatial location; The normalized probability distribution mentioned above is introduced as prior information into the Bayesian inference framework and, together with the borehole hard constraint, drives the interpolation field update of the implicit model, thereby using geophysical data to constrain model inference in areas without borehole coverage.

[0012] Furthermore, the implicit potential energy functional is composed of a weighted sum of three parts: The data fitting term constrains the implicit function to strictly pass through all borehole hard data points and approximates the probability distribution of geophysical soft data in the sense of minimizing residuals; The smoothing term, based on the modulus of the implicit function gradient or the Laplacian operator, controls the geometric smoothness of the generated surface and prevents overfitting. The knowledge violation penalty item is obtained by setting up sampling points in the entire field or a specific area of ​​interest, calculating the potential field deviation of the implicit potential field function at each sampling point that violates the hierarchical rules, morphological rules or topological rules described in the above method, and then summing all the violation deviations by weight. The energy functional achieves a coordinated balance between data fitting accuracy, surface smoothness, and geological logic compliance by adjusting the weight coefficients of the smoothing term and the knowledge violation penalty term.

[0013] Furthermore, the iterative optimization algorithm employs either iterative reweighted least squares or a penalty function method, and the solution process is divided into two stages: In the initial iteration phase, the weight coefficient of the knowledge violation penalty item is set to a small value, so that the algorithm can quickly obtain a geometric approximate solution based on borehole and geophysical data. In subsequent iterations, the weight coefficient of the knowledge violation penalty term is gradually increased, forcing the implicit function to undergo gradual deformation while maintaining the data fitting quality, until the potential field deviation that violates the semantic rules is reduced to below an acceptable threshold, the total energy functional converges, and finally the implicit function solution that satisfies both data constraints and geological logic constraints is output.

[0014] Furthermore, the specific steps for constructing the hierarchical spatial index structure and performing local dynamic updates are as follows: Pre-construct an octree, KD tree, or multi-level volume element index structure for the entire modeling domain to establish a fast retrieval mapping between spatial locations and corresponding implicit equation nodes; When new exploration data is entered into the database, the type of logical conflict between the new data and the existing model is determined based on the engineering geological knowledge graph. If the conflict is determined to be a geometric parameter correction type conflict, the affected spatial radius is determined with the new data coordinates as the center, and the index structure is used to quickly retrieve all nodes that fall within the affected radius. The implicit equations are reconstructed and solved only for the nodes retrieved above, while keeping the potential field parameters of the nodes outside the influence radius unchanged; Through the aforementioned incremental local reconstruction, the three-dimensional geological model can be rapidly corrected without requiring a full recalculation.

[0015] A geological implicit modeling and multi-source data dynamic interpolation system under knowledge graph constraints includes: The data acquisition and knowledge graph management module is used to collect multi-source geological exploration data of the target area and construct or call up engineering geological knowledge graphs; The semantic parsing and mapping module is used to parse the geological semantic rules in the engineering geological knowledge graph, and to transform the sequence relationship rules, geometric morphology rules and structural topology rules into the boundary conditions, anisotropy parameters and penalty functions of the implicit potential field function, respectively. The multi-source data fusion module is used to perform knowledge graph-based confidence assessment on geophysical inversion data and to fuse borehole hard constraints and geophysical soft constraints into a unified probabilistic interpolation field based on Bayesian inference theory. The energy functional construction and solution module is used to construct an implicit potential energy functional that includes data fitting terms, smoothing terms, and knowledge violation penalty terms, and solve the implicit function through an iterative optimization algorithm to extract isosurfaces and generate a three-dimensional geological model. The local dynamic update module is used to maintain a hierarchical spatial index structure covering the entire area. When new exploration data is added to the database, the affected area is quickly located and local incremental reconstruction is performed to realize the real-time dynamic correction of the three-dimensional geological model. The data acquisition and knowledge graph management module, semantic parsing and mapping module, multi-source data fusion module, energy functional construction and solution module, and local dynamic update module work together to execute the method of any one of claims 1 to 8.

[0016] A computer-readable storage medium having a computer program stored thereon, characterized in that, when the computer program is executed by a processor, it implements the method described in any one of the above methods, or drives the system described above to operate.

[0017] Compared with the prior art, the present invention has the following beneficial effects: 1. This invention constructs an engineering geology knowledge graph and establishes a "semantic mathematics" mapping mechanism, transforming geological a priori rules such as stratigraphic sequence, geometric morphology, and tectonic topology into boundary conditions, anisotropic parameters, and penalty functions in the implicit potential field energy functional. The knowledge violation penalty term forces the implicit function to gradually satisfy geological semantic constraints during the iterative solution process, effectively mitigating logical errors such as sequence inversion, morphological distortion, and missing fault movement caused by the "data-only" approach of traditional interpolation algorithms. This ensures that the model generated in sparse borehole areas still conforms to geological genesis, significantly improving model credibility.

[0018] 2. This invention, based on a Bayesian inference framework, defines high-precision borehole data as deterministic hard constraints. Simultaneously, it utilizes geophysical method applicability rules and lithological property mapping rules from a knowledge graph to assess the confidence level of geophysical inversion data, eliminating noisy regions and transforming it into probabilistic soft constraints. This mechanism achieves quantitative fusion of borehole and geophysical data within the same mathematical equation, effectively utilizing geophysical information to constrain the inference range between boreholes. This improves model accuracy while avoiding interference from low-quality data, fully leveraging the value of multi-source exploration data.

[0019] 3. This invention constructs a hierarchical spatial index structure (octree, KD tree, etc.) covering the entire domain. When new exploration data is added to the database, the system quickly determines the conflict type and influence radius based on the knowledge graph. Implicit equations are reconstructed and solved only for the index nodes within the affected area, while external node parameters remain unchanged. This incremental reconstruction strategy reduces model update time from "hours" to "seconds," significantly improving the low efficiency of traditional implicit modeling's full-field recalculation. It enables real-time dynamic evolution of the geological model during construction, providing support for the dynamic updating of the digital twin base plate.

[0020] 4. The implicit functions output by the method and system of this invention can not only extract isosurfaces for visualization, but their mathematical expression and grid units are also naturally compatible with the physical computation grid. Furthermore, the semantic rules in the knowledge graph establish a precise mapping with the spatial units. When engineering conditions change or monitoring data feedback is abnormal, the system can automatically invoke semantic rules to guide local model correction and calculation parameter adjustment, achieving a deep integration of geometric visualization, physical computation, and engineering management. This provides full-lifecycle geological foundation twin support for fields such as water conservancy and hydropower, tunnels, and mining.

[0021] 5. The core algorithm of this application is implemented based on open-source numerical computing libraries (NumPy, SciPy) and graph databases (Neo4j, etc.), which are decoupled from the computing platform and have good cross-platform portability, making it easy to integrate and deploy in various engineering digital platforms. Attached Figure Description

[0022] To more clearly illustrate the technical solutions in the embodiments of this application, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0023] Figure 1 This is a flowchart illustrating a geological implicit modeling and multi-source data dynamic interpolation method according to the present invention. Figure 2 A schematic diagram illustrating the principle of transforming semantic rules (hierarchical order, geometric morphology, and topological construction) into mathematical constraints; Figure 3 This diagram illustrates Bayesian fusion and confidence assessment of multi-source data (borehole hard data and geophysical soft data). Figure 4 This is a schematic diagram illustrating the effect of using a local octree index for dynamic model updates. Figure 5 This is a detailed flowchart of the energy functional construction and iterative solution in the method of the present invention; Figure 6 This is a flowchart of the index retrieval and local reconstruction process with local dynamic updates in the method of this invention; Figure 7 This is a schematic diagram illustrating the module composition and collaborative relationship of a geological implicit modeling and multi-source data dynamic interpolation system according to the present invention. Detailed Implementation

[0024] In the following description, specific details such as particular system architectures and techniques are set forth for illustrative purposes and not for limitation, in order to provide a thorough understanding of the embodiments of this application. However, those skilled in the art will understand that this application may also be implemented in other embodiments without these specific details. In other instances, detailed descriptions of well-known systems, apparatuses, circuits, and methods have been omitted so as not to obscure the description of this application with unnecessary detail.

[0025] This application provides a method and system for implicit geological modeling and dynamic interpolation of multi-source data, which is particularly suitable for engineering scenarios requiring high-precision geological modeling, multi-source data fusion, and dynamic evolution assessment during the design, construction, and operation phases of large-scale water conservancy and hydropower projects (such as high dam foundations, deep-buried tunnels, and large underground powerhouses), as well as transportation tunnels and deep mines. During the construction and operation of large-scale water conservancy and hydropower projects, the spatial morphology and physical property distribution of complex geological structures (such as faults, weak interlayers, dikes, and fracture zones) directly determine the stability of the surrounding rock and the safety of engineering construction. The deformation and failure of the geological substrate are essentially the result of the coupling effect of geometric morphology, semantic knowledge, and geomechanical parameters. There is a strong correlation between the geometric boundary morphology of the geological body, stratigraphic sequence and structural topological relationships (semantic knowledge), and rock mass mechanical parameters. Simplification or fragmentation of any link will significantly affect the reliability of the numerical analysis results. However, traditional 3D implicit geological modeling methods suffer from three major flaws: First, relying solely on geometric distance interpolation often generates geologically inconsistencies (such as sequence inversions and missing fault movements) in sparse borehole areas, making them mathematically driven models that "understand mathematics but not geology." Second, it's difficult to scientifically integrate borehole hard data and geophysical soft data within the same mathematical equation, leaving a large amount of geophysical data idle and lacking effective constraints on the model's predictive range. Third, implicit modeling involves solving large systems of linear equations, requiring a complete recalculation after new data is revealed during construction, resulting in significant model update delays. This makes it difficult to dynamically update the geological model with the construction progress and achieve dynamic evolution during construction and operation. These shortcomings lead to significant discrepancies between the geological model's predicted mechanical response and actual monitoring results, making it difficult to effectively identify geological risk areas induced by construction excavation and affecting the scientific formulation of engineering support schemes and safety monitoring measures.

[0026] This invention constructs an engineering geological knowledge graph and establishes a "semantic mathematics" mapping mechanism to transform geological rules such as sequence stratigraphy, morphology, and structure into implicit potential field function boundary conditions, anisotropic parameters, and penalty functions. Based on a Bayesian inference framework, borehole data is defined as deterministic hard constraints, while geophysical inversion data is transformed into probabilistic soft constraints after knowledge graph confidence assessment, achieving multi-source data fusion. An energy functional containing data fitting terms, smoothing terms, and knowledge violation penalty terms is constructed, and the implicit functional expression of geological interfaces is solved using iterative reweighted least squares or penalty function methods. Furthermore, a hierarchical spatial index structure (octree, KD tree, etc.) covering the entire domain is constructed. When new exploration data is added to the database, the conflict type is determined based on the knowledge graph, and the influence radius is determined. Only the index nodes in the affected area are locally incrementally reconstructed, achieving second-level dynamic updates of the model. The aforementioned technical approach fundamentally solves the three major problems in existing implicit modeling: lack of geological logic, difficulty in multi-source data fusion, and lag in dynamic updates. It provides high-precision, logically strong, and fast-response digital twin support for the geological base of water conservancy and hydropower projects throughout their entire life cycle for intelligent construction and safety decision-making.

[0027] In this embodiment, the engineering geology knowledge graph is constructed using the Python programming language combined with open-source numerical computing libraries (such as NumPy and SciPy) and graph databases (such as Neo4j or JanusGraph). Radial basis functions (RBF) or Kriging interpolation are used as basis functions for the implicit potential field function. The energy functional is solved by iterative reweighted least squares (IRLS) method, and local dynamic updates are achieved based on an octree index structure. This forms a complete technical closed loop from geological semantic parsing, Bayesian fusion of multi-source data, implicit potential field solution to local incremental updates. The technical solution of this application will be described in detail below with reference to specific embodiments.

[0028] Example 1 Embodiment 1 of this application provides a method for implicit geological modeling and dynamic interpolation of multi-source data, such as... Figure 1 As shown, the method includes steps S101 to S105. Example 1 constructs an engineering geology knowledge graph and establishes a semantic-mathematical mapping mechanism, transforming sequence stratigraphy, geometric morphology, and structural topology rules into implicit potential field function boundary conditions, anisotropic parameters, and penalty functions; it integrates borehole hard constraints and geophysical soft constraints based on Bayesian inference; it constructs and solves an energy functional containing data fitting terms, smoothing terms, and knowledge violation penalty terms; and finally, it uses a hierarchical spatial index to achieve local dynamic updates. This method solves three major technical problems in existing implicit modeling: missing geological logic, difficulty in multi-source data fusion, and lagging dynamic updates. The specific implementation is as follows: S101. Obtain multi-source geological exploration data for the target area and construct or call up an engineering geological knowledge graph: In application, acquiring multi-source geological exploration data of the target area is fundamental to building a three-dimensional geological model. This process includes collecting borehole data, such as hard data like lithology, stratigraphy, and rock mechanical parameters, as well as geophysical data such as resistivity and wave velocity inversion results. The engineering geological knowledge graph is a structured geological information database containing geological entities such as strata, faults, and rock masses; geological attributes such as age, occurrence, and hardness; and semantic relationships between entities such as overlying relationships and cutting relationships. In the dam site modeling, hard data from 20 boreholes (ZK1 to ZK20) were collected, along with controlled-source audio-frequency magnetotelluric (CSAMT) geophysical profile data. The constructed engineering geological knowledge graph records geological entities such as Quaternary overburden, granite dikes, and sedimentary bedrock, as well as their sequence relationships, genetic characteristics, and tectonic cutting relationships.

[0029] As an example, this embodiment uses the Resource Description Framework (RDF) standard to construct a triple data structure for "entity-relationship entities". The core stored content of the knowledge graph is listed below: (1) Triple with ID 001: Subject is Geo:Q4 (fourth series), predicate is rdf:type, object is Geo:Sedimentary_Layer, semantic parsing is to define Q4 as a sedimentary layer entity.

[0030] (2) The triple with ID 002: the subject is Geo:T (sedimentary rock), the predicate is rdf:type, the object is Geo:Bedrock, and the semantic parsing is to define T as a bedrock entity.

[0031] (3) Triple with ID 003: Subject is Geo:Q4, predicate is rel:overlies, object is Geo:T, semantic parsing is geological rule A:Q4 is above T (source of inequality constraint).

[0032] (4) Triple with ID 004: Subject is Geo:γ (dipole), predicate is rdf:type, object is Geo:Intrusive_Vein, semantic parsing is to define γ as an intrusive dipole.

[0033] (5) Triple with ID 005: subject is Geo:γ, predicate is prop:shape_type, object is Shape:Tabular, semantic parsing is triggering anisotropic kernel function selection.

[0034] (6) Triple with ID 006: subject is Geo:γ, predicate is prop:main_axis_ratio, object is 10:10:1 (main:secondary:method), semantic parsing is geological rule B: defines the anisotropy ratio of the covariance matrix.

[0035] (7) Triple with ID 007: Subject is Geo:F1 (fault), predicate is rel:cuts (cut), object is Geo:T, semantic parsing is triggering fault cutting relationship processing.

[0036] (8) Triple with ID 008: Subject is Geo:F1, predicate is rel:drags (traction), object is Geo:T, semantic parsing is triggering local orientation rotation matrix calculation.

[0037] (9) Triple with ID 009: subject is Geo:F1_Drag, predicate is prop:rotation_angle, object is 15_degrees, semantic parsing is geological rule C: set rotation parameters.

[0038] The system reads the above triples using the SPARQL query language and automatically activates the corresponding mathematical constraints in the energy functional.

[0039] S102. Analyze the geological semantic rules in the knowledge graph and establish a "semantic-mathematical" mapping mechanism, such as... Figure 2 As shown, the system transforms the hierarchical sequence, geometric shape, and construction topology rules into inequality constraints, anisotropy matrices, and step / rotation matrices, respectively: In application, parsing the geological semantic rules in the engineering geology knowledge graph is key to converting geological knowledge into mathematical constraints. Stratigraphic sequence rules describe the vertical relationships between different geological bodies in space, such as Quaternary overlying granite and sedimentary rocks. Geometric morphology rules describe the genetic characteristics and occurrence of geological bodies, such as dikes extending in a plate-like pattern, dipping at NE30° with a dip angle of 75°. Tectonic topology rules describe the cutting and traction relationships between tectonic bodies and other geological bodies. Establishing a mapping relationship between semantic rules and mathematical constraints enables the computer to understand and execute these geological logics. In the dam site implementation example, the system adopts the resource description framework standard and constructs a triplet data structure of entity-relationship-entity, allowing geological rules such as Q4 overlying T, dikes being plate-like, and fault F1 cutting through T to be recognized by the computer system.

[0040] In application, sequence relation rules are transformed into boundary conditions for implicit potential field functions. When the knowledge graph defines Quaternary strata as overlying sedimentary bedrock, the system constructs inequality constraints in the normal direction of the stratigraphic contact surface, forcing the potential field value of the Quaternary strata to be greater than that of the sedimentary bedrock at any spatial point. This constraint prevents stratigraphic inversion caused by mathematical smoothing in sparsely drilled areas. Geometric morphology rules are transformed into anisotropic parameters, constructing covariance matrices or weight matrices based on sedimentary or tectonic origins. For stratiform sedimentary geological bodies, the matrix ensures that the range of the interpolation kernel function along the strike-dip direction is significantly greater than its range along the normal direction, allowing lenses or thin layers to exhibit correct geometric morphology. For vein-intrusive geological bodies, the matrix makes the kernel function flattened ellipsoidal, consistent with the dike attitude. Tectonic topology rules are transformed into penalty functions, introducing step functions or drift terms related to the fault distance on both sides of the fault in the implicit equations to express the displacement discontinuity of the strata on both sides of the fault.

[0041] Regarding the mathematical expression of fault cutting relationships, this embodiment provides complete definitions of the drift term and the trend rotation matrix: (1) Fault plane sign distance function Let the equation of the fault plane be Φ(x)=0, where Φ(x) is the sign distance function, satisfying that |Φ(x)| is equal to the distance from point x to the fault plane, Φ(x)>0 represents the hanging wall, and Φ(x)<0 represents the footwall.

[0042] (2) Step function Introducing a step function H(Φ) to express the displacement discontinuity on both sides of the fault: Standard step form: When Φ>0, H(Φ)=1; When Φ ≤ 0, H(Φ) = 0; Smooth approximation (to avoid numerical instability): Where σ is a smooth transition parameter that controls the gradual range of displacement near the fault.

[0043] (3) Drift items with rotation The drift term, which includes displacement and orientation rotation, is defined as follows: Where d(x) is the direction-dependent displacement vector, its expression is: Meaning of symbols in the formula: v: The basic displacement vector of the fault, which is determined by the motion properties and displacement of the fault; R(θ): Rotation matrix around the fault strike, used to simulate the rotation of the hanging wall relative to the footwall. L: Displacement attenuation length, which controls the rate at which the displacement gradually decreases from its maximum value at the fault plane to both sides; θ: Rotation angle, determined based on the traction structure parameters recorded in the knowledge graph.

[0044] (4) The specific form of the rotation matrix R(θ) R(θ) = [1, 0, 0; 0, cosθ, -sinθ; 0, sinθ, cosθ] Right now: First line: [1, 0, 0]; The second line: [0, cosθ, -sinθ]; The third line: [0, sinθ, cosθ]; (5) Parameter values ​​in this embodiment The rotation angle θ = 15° (derived from the prop:rotation_angle attribute of Geo:F1_Drag in the knowledge graph). Displacement attenuation length L = 20 m (set according to the range of the fault influence zone); The smooth transition parameter σ = 2m; This rotation matrix acts on the dip and normal directions in the local coordinate system of the fault plane, causing the attitude of the strata on the hanging wall to rotate by 15° relative to the footwall, thereby simulating the traction deformation associated with fault activity. Based on the motion characteristics of the fault, a trend rotation matrix is ​​added to the drift term to simulate the attitude rotation and traction deformation of the hanging wall relative to the footwall.

[0045] Regarding the construction of the anisotropic matrix described above, this embodiment provides its explicit mathematical definition and weighted distance calculation method: (1) Construction of rotation matrix R: The geological body has a strike of α and a dip angle of β. In this embodiment, the dike has a strike of α = 30° (NE30°) and a dip angle of β = 75°. The global coordinate system (x, y, z) is transformed to the local coordinate system (x′, y′, z′), where the x′ axis is along the strike direction, the y′ axis is along the dip direction, and the z′ axis is along the normal direction. The rotation matrix R is defined as: Substituting α=30° and β=75° yields the specific numerical matrix.

[0046] (2) Construction of the anisotropic variable range matrix Λ set up , , These represent the range (affected area) along the strike, dip, and normal directions. Based on the "plate-like" geometry of the dike recorded in the knowledge graph, the range ratio is set as follows: : : =10:10:1. The range matrix Λ is defined as: (3) Construction of the anisotropic covariance matrix M: Combining the rotation matrix and the range matrix yields the covariance matrix: This matrix is ​​a third-order symmetric positive definite matrix used to define anisotropic distances.

[0047] (4) Calculation of anisotropic weighted distance: For any two points in space, u = ( , , ) and v=( , , The anisotropic distance is defined as: Equivalently, calculated in the local coordinate system: in .

[0048] (5) Anisotropic radial basis functions Gaussian radial basis functions are used as the interpolation kernel function for the implicit potential field: in For global smoothing parameters, the following values ​​can be taken: =max( , , This kernel function ensures that the data influence along the strike and dip directions of the vein is much greater than that along the normal direction, thus guaranteeing that the generated vein model presents a flattened ellipsoidal shape, consistent with the true geological morphology of plate-like extension.

[0049] S103. Construct a multi-source data fusion field based on Bayesian inference, such as Figure 3 As shown, the system performs confidence assessment on geophysical anomaly areas, eliminates low-confidence noise, and generates probabilistic soft constraints: In the application, Bayesian inference theory is used to fuse multi-source geological exploration data. Borehole exploration data, derived from deep excavation, is highly accurate and reliable, and is defined as a deterministic hard constraint, meaning the probability of a geological attribute existing at the borehole coordinate point is 100%. Geophysical inversion data, obtained through inversion of stratigraphic structures using physical fields such as electromagnetics and seismic data, contains multiple solutions and noise, requiring confidence assessment before being converted into probabilistic soft constraints. In the dam site example, low resistivity anomalies may correspond to fault fracture zones or karst areas. The system's knowledge graph query revealed that the bedrock in this area is dense limestone. The location of the fracture zone revealed by borehole ZK5 is used to verify the geophysical anomaly's consistency, and the system calculates the overall confidence level. in The applicability coefficient of the geophysical method. The property-lithology mapping matching degree, where =0.8, based on the applicability rating of the CSAMT method in tight limestone areas in the knowledge graph; =0.9, calculated based on the spatial agreement (85%) between the low resistivity anomaly and the fractured zone revealed by borehole ZK5. Transforming geophysical anomalies into probability distributions of geological attributes effectively constrains model predictions for areas without borehole coverage with geophysical data.

[0050] The specific form of the Bayesian update is as follows: Let P(G) be the prior probability of the geological attribute determined based on borehole hard data (1 at the borehole point, and calculated by the interpolation kernel function at other locations), and P(D|G) be the likelihood function of the geophysical anomaly data under a given geological attribute (given by the property-lithology mapping rule). Then, the posterior probability after geophysical observation is P(G|D) ∝ P(D|G)·P(G). This invention uses the normalized confidence coefficient Cscore of high-confidence geophysical anomaly areas to approximate the contribution of the likelihood function, and achieves quantitative fusion of prior and geophysical data through a multiplicative update rule.

[0051] S104. Construct and solve the implicit potential energy functional, generate a three-dimensional geological model, and implement local dynamic updates, such as... Figure 5 As shown, a two-stage IRLS iterative optimization is adopted, with the initial... It is relatively small, and gradually increases in size: In the application, an implicit potential energy functional is constructed, comprising a data fitting term, a smoothing term, and a knowledge violation penalty term. The data fitting term constrains the implicit function to strictly pass through all borehole hard data points and approximates the probability distribution of geophysical soft data in the sense of minimizing residuals. For borehole hard data points, when constructing the design matrix A, the corresponding rows are assigned weight coefficients that are much higher than those for soft data (in this embodiment, we take...). = 1×10 6), to approximately achieve the equality constraint f( within the range of real number calculation precision. ) = While maintaining the uniformity of the overall optimization problem form, it facilitates the use of standard least squares solvers. The smoothing term, based on the modulus or Laplacian operator of the implicit function gradient, controls the geometric smoothness of the generated surface, preventing overfitting. The knowledge violation penalty term is obtained by calculating the potential field deviation at each sampling point (where the implicit potential field function violates stratigraphic, morphological, or topological rules) by deploying sampling points across the entire field or within a specific area of ​​interest, and then summing all the violation deviations using a weighted average. By adjusting the weight coefficients of the smoothing term and the knowledge violation penalty term, a coordinated balance is achieved between data fitting accuracy, surface smoothness, and geological logical compliance. In the dam site example, the overall objective function is expressed as: in: This is a data fitting term used to constrain the implicit function to pass through hard data points and approximate the probability distribution of soft data; For the smoothing term, the gradient-based modulus or Laplacian operator controls the surface smoothness; For knowledge violation penalty items, the weighted sum of potential field deviations that violate the semantic rules of step S102 is calculated by sampling the entire field or a specific region. , These are the weighting coefficients.

[0052] In the application, an iterative optimization algorithm is used to solve the energy functional. In the initial iteration phase, the weight coefficient of the knowledge violation penalty term is set to a small value, allowing the algorithm to quickly obtain a geometric approximation solution based on borehole and geophysical data. At this stage, the model may contain areas that violate geological rules, such as dikes penetrating Quaternary overburden. In subsequent iterations, the weight coefficient of the knowledge violation penalty term is gradually increased, forcing the implicit function to undergo asymptotic deformation while maintaining the data fitting quality, until the potential field deviation that violates semantic rules is reduced to below an acceptable threshold, and the total energy functional converges.

[0053] Penalties for knowledge violations For specific calculations, this embodiment uses the sequence rule of stratigraphy as an example to give its mathematical definition. Let the implicit potential field value of stratigraphy A at the spatial sampling point xi (i=1,2,...,M) be... The implicit potential field value of stratum B is According to the stratigraphic sequence rule defined in the knowledge graph, "stratum A overlies stratum B", the following conditions should be met: > Define the penalty for a single violation. for: Summing the violation penalties for all M sampling points yields the total knowledge violation penalty term: In the total energy functional, this term is multiplied by the weighting coefficient. Later, we participated in the optimization: This embodiment uses the Iterative Reweighted Least Squares (IRLS) method to solve the above energy functional. The core idea is to transform the original problem into a series of weighted least squares subproblems, and then approximate the optimal solution of the original objective function by iteratively updating the weight matrix. The specific iterative steps are as follows: (1) Initialization: Set the iteration counter k=0, the initial coefficient vector x(0) (which can be the zero vector or the least squares solution based on the borehole hard data), and the initial knowledge violation penalty term weight coefficient. =0, weighting factor γ=1.5, maximum weighting factor =100, convergence threshold ε= Smoothing term weight coefficient Based on a pre-set data noise level (in this embodiment, it is taken as...). =0.1).

[0054] (2) Constructing a weighted least squares problem: In the k-th iteration, based on the current solution Calculate the weight of each sampling point .

[0055] For soft data points in the data fitting term, the weights are taken as confidence levels. = (Given by the confidence assessment in step S103).

[0056] For the sampling points in the knowledge violation penalty item, the weights are approximated by the first derivative of the penalty function: Where δ= To prevent small positive numbers from being divided by zero.

[0057] (3) Solve the weighted least squares subproblem: Solve the following system of linear equations to obtain a new coefficient vector. : in: A is the design matrix, which consists of the values ​​of each radial basis function at the sampling points; This is the weight diagonal matrix for the current iteration step; L is the discretization matrix of the smoothing term (such as the Laplace matrix). The Hessian matrix of the knowledge violation penalty term is composed of the second derivatives of the violation sampling points; b is the observation vector (normalized values ​​of borehole hard data and geophysical soft data).

[0058] (4) Update the weight coefficients of knowledge violation penalty items: adopt an exponential increasing strategy: (5) Convergence criterion: Calculate the change in solution between two consecutive iterations: If Δ < ε and ≥ If it converges, then the output will be... Otherwise, let k = k + 1 and return to step (2) to continue the iteration.

[0059] In the above iteration process, steps (2)-(3) transform the original nonlinear nonconvex optimization problem into a series of linear least squares subproblems, and step (4) gradually increases the... This forces the solution to converge in a direction that satisfies geological semantic rules. The algorithm achieves a balance between data fitting and geological logic compliance while ensuring computational stability.

[0060] In the dam site area embodiment, the iterative reweighted least squares method is used for solution. In the initial iteration, only borehole and geophysical data are fitted. In the intermediate iteration, sampling detects the violation point where the dike penetrates the Quaternary system, and the penalty term increases sharply. The algorithm automatically reduces the potential field value of the dike. Finally, the model is obtained where the dike is truncated at the bottom boundary of the Quaternary system, has a plate-like shape, and the strata on both sides of the fault are clearly displaced.

[0061] like Figure 4 and Figure 6 As shown, octree indexing is used to quickly locate nodes within the affected region of interest (ROI), with only local reconstruction, achieving updates within seconds. To ensure the continuity of the updated model at the ROI boundaries, this invention employs the following mechanism: The radius of influence is determined by the distance at which the radial basis function kernel value decays to less than 1% of the peak value. Within this range, the contribution of the basis function is close to zero, ensuring that the parameters of the nodes outside the boundary remain unchanged and do not introduce discontinuity breaks. For nodes near the ROI boundary (buffer zone less than 5m from the ROI boundary), a linearly gradient weight is used to mix the old and new solution results to ensure that the gradient of the potential function is continuous at the boundary.

[0062] In application, isosurfaces are extracted based on the solved implicit functions to generate a 3D geological model. Classical algorithms such as the moving cubes are used to extract the isosurfaces of the implicit functions, i.e., finding all points in space that satisfy the implicit equation f(x,y,z) equal to a specific potential field value, and connecting them to form a geological interface in the form of a triangular mesh. A hierarchical spatial index structure is constructed, using octrees, KD trees, or multi-level voxel indexes to establish a fast retrieval mapping between spatial locations and corresponding implicit equation nodes. When newly added exploration data arrives, the system first determines the type of logical conflict between the new data and the existing model based on the engineering geology knowledge graph. If it is determined to be a geometric parameter correction conflict, the affected spatial radius is determined with the new data coordinates as the center. The index structure is used to quickly retrieve all nodes falling within this radius, and the implicit equations are reconstructed and solved only for these nodes, keeping the potential field parameters of nodes outside the radius unchanged. Through this incremental local reconstruction, the 3D geological model can be quickly corrected without requiring a full-field recalculation. In the dam site construction example, when the tunnel boring machine (TBM) reached chainage K0+500, the actual thickness of the exposed rock vein was 2 meters less than the predicted thickness. The system quickly located the volume element nodes within a 30-meter radius of this location using an octree index, and only resolved the implicit equations for the internal nodes, completing the local model correction in 0.5 seconds. Experimental tests show that in a typical dam site modeling scenario containing 100,000 nodes, the traditional full-site recalculation scheme takes an average of about 4.2 hours, while the local incremental reconstruction scheme of this invention (ROI radius 30m, containing about 200 affected nodes) takes an average of about 0.5 seconds, improving efficiency by about 3×10⁻⁶. 4 The stratigraphic sequence violation rate was reduced from 18.3% with traditional RBF interpolation to less than 0.2% with this method, verifying the effectiveness of knowledge constraints.

[0063] Specifically, the entire modeling process forms a closed-loop dynamic evolution system. Based on the acquisition of multi-source geological data, the system transforms the prior experience and rules of geological experts into mathematical constraints through a knowledge graph, making implicit modeling no longer purely data-driven, but an organic integration of knowledge and data. In the multi-source data fusion stage, the Bayesian inference framework enables data of varying precision to be scientifically integrated according to their respective credibility levels. Hard data has strong constraints, while soft data, filtered by the knowledge graph, becomes an effective supplementary constraint. In the energy functional solution stage, by progressively increasing the weight of knowledge violation penalty terms, the optimization algorithm gradually adjusts towards satisfying geological logic while pursuing data fitting accuracy, ultimately obtaining a 3D model that conforms to both observational data and geological laws. In the dynamic update stage, the hierarchical spatial index structure allows the introduction of new data to avoid full-field recalculation, requiring only incremental reconstruction in locally affected areas, significantly improving update efficiency. This working principle ensures the logical consistency of the model, the scientific nature of data fusion, and the high efficiency of updates.

[0064] The stratigraphic sequence constraints can be implemented by constructing inequality constraints along the normal direction of the stratigraphic contact surface, forcing the potential field value of the overlying strata to always be greater than that of the underlying strata. Alternatively, a secondary penalty term can be applied when calculating the negative value of stratigraphic thickness at all sampling points across the entire field, ensuring that non-geological stratigraphic inversions do not occur in sparse data regions. The anisotropy parameter construction can utilize a covariance matrix based on the strike, dip, and dip angle of the geological body. By changing the eigenvectors and eigenvalues ​​of the matrix, the extension direction and range ratio of the interpolation kernel function can be controlled. Alternatively, predefined standard anisotropy templates based on the genetic type of the geological body (sedimentary, intrusive, metamorphic, etc.) in the knowledge graph can be used, with the system automatically selecting the best-matching template for rapid configuration. The iterative solution of the energy functional can employ iterative reweighted least squares, updating the objective function for each iteration by calculating the L2 norm of the weighted residuals. A penalty function method can also be used, transforming the constrained optimization problem into a sequential solution of an unconstrained problem by introducing an external penalty function. The retrieval of the index structure can be implemented using an octree, which recursively divides the space into eight parts, with each leaf node corresponding to a voxel. During the query, the query proceeds from the root down along the tree structure to find the voxel where the target point is located and its adjacent voxels. Alternatively, a KD tree can be used, which constructs a binary tree structure by alternating the partitioning of the space in different dimensions, enabling fast range queries and nearest neighbor searches.

[0065] In terms of technical advantages, the method framework established in this claim solves the problem of traditional implicit modeling, which purely pursues mathematical smoothness while neglecting geological logic. By transforming geological semantic rules into constraint and penalty terms in energy functionals, the generated 3D model can conform to geological genesis laws even in sparse data regions, significantly improving the reliability of predictions. The Bayesian fusion framework for multi-source data scientifically handles the heterogeneity of hard and soft data, making full use of the constraints of geophysical data on borehole-free areas and eliminating the interference of geophysical noise through the filtering mechanism of knowledge graphs, thus improving the accuracy of the overall model. The combination of hierarchical spatial index structure and local incremental reconstruction algorithm reduces the model update time from hours to seconds, supporting the real-time correction needs during construction and truly realizing the dynamic evolution of the geological model as construction progresses, meeting the high-efficiency adaptability requirements of digital twin base plates. In addition, the core algorithm of this method is independently controllable and does not rely on foreign commercial modeling software, which is conducive to protecting the security of national basic geological data.

[0066] As an optional implementation, the specific way to transform the hierarchical relationship rule into boundary conditions is as follows: In applications, stratigraphic sequence rules define the temporal and spatial relationships between different geological bodies. In geology, according to the stratigraphic overlay law, the overlying strata are deposited later than the underlying strata; this strict vertical sequence relationship is the basis for stratigraphic division and correlation. In knowledge graphs, this relationship is explicitly encoded as rel:overlies, for example, Q4 overlying T indicates that Quaternary sediments accumulated on sedimentary bedrock. The potential field value of the implicit potential field function usually represents the relative age of the strata, with a larger potential field value corresponding to a newer stratum. Therefore, the stratigraphic interface is where the potential field values ​​of two adjacent strata change most drastically; along the normal direction of the interface, the potential field gradient direction of the overlying strata should point towards the more recent strata.

[0067] In application, the specific way to transform the stratigraphic sequence rules into boundary conditions is to construct inequality constraints in the normal direction of the stratigraphic contact surface. If the knowledge graph defines stratigraphy A as overlying stratigraphy B, then the system adds constraint formulas to the implicit modeling constraints, requiring the implicit potential field value of stratigraphy A at any spatial point. (x) is greater than the implicit potential field value of stratum B. (x), that is (x) - (x)>0. This constraint aims to prevent non-geological stratigraphic inversions during the modeling process. In areas with sparse borehole data, simply interpolating based on geometric distance may produce physically unreasonable results, such as the potential field value of the upper Quaternary overburden being smaller than that of the lower bedrock, leading to morphological stratigraphic inversion. Inequality constraints, by enforcing constraints, ensure that the optimization algorithm, while satisfying data fitting, guarantees the correct sequence of potential field values ​​between the upper and lower strata.

[0068] In the application, the specific scenario of geological modeling in the dam site area involves three main strata: Quaternary, granite dikes, and sedimentary bedrock. The knowledge graph clearly records the stratigraphic sequence, with Q4 overlying γ and T, meaning the Quaternary is at the top and the bedrock at the bottom. The system identifies constraint points at all stratigraphic contact surfaces traversed by boreholes; for example, borehole ZK3 crosses the contact surface between Q4 and T at a depth of 25 meters. Inequality constraints are set in the normal direction near this contact surface, requiring that the potential field value of Q4 is greater than the potential field value of T throughout the entire modeling area. When a sampling point violates this constraint during iterative optimization, such as a location where the potential field value of Q4 is less than or equal to the potential field value of T, a penalty term is activated at that point, generating a large energy increment that forces the optimization algorithm to adjust relevant parameters to eliminate the violation. After multiple iterations, the final generated model shows a clear interface between Q4 and T, ensuring the geological rationality of the stratigraphic sequence.

[0069] Specifically, the implementation of stratigraphic sequence constraints in implicit modeling follows these principles. First, the system retrieves all overriding relationships from the knowledge graph. For each ordered stratigraphic pair (A, B), where A overrides B, the system adds a constraint term to the implicit potential energy functional. This constraint term is not defined on a two-dimensional surface at the contact point, but rather across the entire three-dimensional modeling region; any spatial location must satisfy the sequence constraint of the potential field values. During the iterative solution process, the system periodically samples and detects, randomly or uniformly scattering sampling points in space, and calculating the magnitude relationship of the potential field values ​​of different stratigraphic layers at each sampling point. If a sampling point is found to be... ≤ This indicates a violation, and the potential field deviation at that sampling point is defined as max(0, - ) 2 The sum of the deviations of all non-compliant sampling points is multiplied by a weighting factor. This is added to the total energy functional. As the weight coefficients gradually increase, the energy contribution of the violation term gradually increases, forcing the optimization algorithm to adjust the parameters of the implicit function so that at each sampling point... > This continues until all sampling points satisfy the constraints.

[0070] Among them, the definition of inequality constraints can adopt strict inequalities. (x) - If (x)>0, and the potential difference is required to always be positive, a buffered inequality can also be used. (x) - (x)>ε, where ε is the minimum separation threshold greater than 0, to maintain sufficient boundary clarity even when potential field values ​​are close. The constraint points can be set densely at all borehole penetration points on the known formation contact surface to ensure constraint strength in key areas, or uniformly sampled along the entire contact surface curve to create coherent constraint conditions even in borehole-free areas. The penalty function can be a quadratic form max(0, - ) 2 This makes the penalties stronger for larger violations, and higher-order forms such as the quartic form max(0, - ) 4 More severe penalties will be imposed for serious violations.

[0071] In terms of technical advantages, by applying stratigraphic sequence constraints throughout the entire modeling domain, this method ensures that the generated 3D model conforms to the true vertical relationship of strata at any location, completely eliminating non-geologically derived stratigraphic inversion phenomena. This constraint is a hard constraint implemented at the mathematical level, which the optimization algorithm must adhere to, rather than requiring post-hoc verification or correction. Compared to the stratigraphic inversion problem that easily occurs in sparsely drilled areas with traditional implicit modeling, this method ensures the inevitability of geological logic through knowledge-driven constraints, significantly improving the geological rationality of the prediction results. For engineering applications, this means that decisions made based on this model have greater credibility and can reduce engineering risks caused by model errors.

[0072] As an optional embodiment, the specific method for converting geometric morphological rules into anisotropic parameters is as follows: In application, geometric morphological rules refer to the spatial morphological characteristics of a geological body based on its genetic type and occurrence. The genetic types of strata include sedimentary, intrusive, and metamorphic origins, and geological bodies exhibit different spatial extension patterns under different genetic backgrounds. Sedimentary strata are often controlled by the sedimentary paleogeographic environment and gravity, exhibiting layered or lenticular forms, mainly extending along the direction of ancient sedimentary flows, with relatively uniform thickness in the strike direction and thinner in the normal direction. Intrusive dikes intrude along the weak surfaces of magma ascents, often exhibiting vein-like or network-like forms, mainly extending along the strike and dip of the vein, while the thickness perpendicular to the vein surface is usually smaller. Occurrence attributes include three parameters: strike, dip, and dip angle, which fully describe the spatial directional characteristics of a geological body.

[0073] In application, the specific method for transforming geometric morphology rules into anisotropic parameters is to construct an anisotropic covariance matrix or weight matrix based on the genetic type and occurrence attributes of geological bodies recorded in the engineering geological knowledge graph. This matrix is ​​a third-order symmetric positive definite matrix, whose eigenvectors define three orthogonal principal directions, and eigenvalues ​​define the range ratio in each direction. Range refers to the influence range of the interpolation kernel function; a larger range indicates a larger influence range for a data point in that direction. This matrix is ​​used to control the range ratio and rotation angle of the interpolation kernel function in different spatial directions, ensuring that the generated geological body morphology conforms to its sedimentary or tectonic genetic background.

[0074] In applications, for sedimentary geological bodies extending in layers, the matrix ensures that the range of the interpolation kernel function along the strike-dip direction is significantly greater than the range along the normal direction. In the embodiment of dam site modeling, the attitude attributes of the granite dike are recorded in the knowledge graph as strike NE30°, dip SE, and dip angle 75°. The anisotropy matrix established by the system based on these parameters sets the ratio of the three principal axes to 10:10:1. The first two ratios of 10 indicate that the range along both the strike and dip directions of the dike is 10 units, and the third ratio of 1 indicates that the range perpendicular to the dike surface (i.e., the normal direction of the dike) is 1 unit. As a result of this setting, when traversing any point in space, the influence of borehole data 100 meters away from that point along the dike's extension direction is comparable to the influence of borehole data only 10 meters away from that point along the normal direction. Therefore, even in areas far from the normal of the vein, information along the vein's extension direction can be effectively transmitted, resulting in a flattened ellipsoidal shape with the major axis along the vein's direction and the minor axis along the vein's normal, perfectly matching the true occurrence characteristics of the vein.

[0075] In application, for intrusive geological bodies occurring in vein-like patterns, the matrix makes the interpolation kernel function flattened into an ellipsoid, consistent with the dike's occurrence. This setting is particularly suitable for modeling intrusive rock masses such as granite dikes and diabase dikes. In the dam site area, when the thickness of the dike exposed during TBM tunnel construction is 3 meters, the system can accurately infer the thickness and extent of the dike at other locations based on the previously established anisotropy matrix, without excessively deforming the dike thickness due to the influence of a single borehole point.

[0076] Specifically, the complete process of transforming geometric morphology rules into anisotropic parameters is as follows. First, the system queries the knowledge graph for the genetic type and attitude attributes of geological bodies. For layered sedimentary strata, the system extracts three attitude parameters: strike, dip, and dip angle, and establishes a coordinate transformation matrix such that the first and second axes in the transformed coordinate system are along the strike-dip direction of the strata, and the third axis is along the normal direction of the strata. Then, an anisotropic matrix is ​​constructed in this coordinate system, setting larger eigenvalues ​​(corresponding to large ranges) in the strike-dip direction and smaller eigenvalues ​​(corresponding to small ranges) in the normal direction. For vein-like intrusive terranes, the system uses a similar method, but the configuration of eigenvalues ​​is more extreme, making the range perpendicular to the vein surface much smaller than the range along the vein surface. Finally, the system transforms this matrix back to the global coordinate system through coordinate transformation, obtaining the final anisotropic covariance matrix used for implicit modeling. In the interpolation kernel function of implicit modeling, this matrix changes the weighted Euclidean distance calculation method of the radial basis function, making the influence of data points on their surrounding space exhibit anisotropic distribution, rather than the traditional isotropic spherical distribution.

[0077] The eigenvalue ratios of the anisotropy matrix can be configured based on the actual ratio of the thickness and extension length of the geological body, ensuring that the range ratio matches the true proportion of the geometric shape. Alternatively, they can be configured based on predefined genetic type standard templates in the knowledge graph; for example, the default setting for sedimentary strata is that the strike range is 10 times the normal range. The rotation of the matrix in the coordinate system can be determined by directly calculating the rotation angle based on the three attitude parameters to identify the three eigenvectors of the matrix. Alternatively, the rotation angle can be automatically optimized using geometric constraints in the knowledge graph to ensure that the matrix orientation accurately corresponds to the actual attitude of the geological body. The interpolation kernel function can be implemented using a Gaussian radial basis function combined with Mahalanobis distance calculation, causing the kernel value to decrease in a Gaussian distribution with anisotropic distance. Thin-plate spline functions combined with a weighted covariance matrix can also be used to achieve smoother surface fitting characteristics.

[0078] In terms of technical advantages, by configuring anisotropic parameters, implicit modeling can automatically adjust the directionality and strength of interpolation based on the genetic characteristics and occurrence attributes of geological bodies, ensuring that the generated geological bodies naturally conform to geological principles without requiring manual post-processing corrections. For layered strata, the model can accurately reproduce the bedding characteristics of the strata, avoiding the erroneous expansion of stratum thickness due to sparse borehole data in the vertical direction. For vein-like rock masses, the model can constrain the width of the veins, preventing excessive widening of the veins due to abundant data along the vein direction. This method is particularly suitable for complex terrain areas with unevenly distributed engineering geological exploration data; through knowledge-driven anisotropic settings, it can significantly improve the prediction accuracy in sparse data areas.

[0079] As an optional implementation, the specific way to transform the construction topology rules into penalty functions is as follows: In application, structural topology rules refer to the spatial relationships between faults, folds, and other geological structures, particularly cutting and traction relationships. A fault is a structure in which strata break and undergo relative displacement along the break surface. Its most prominent characteristic is that adjacent strata on both sides of the fault plane are displaced, disrupting the continuity of the strata. According to Stino's law, under normal circumstances, strata are horizontally deposited, but fault activity causes strata at different locations to shift. Simultaneously, fault movement is often accompanied by strata rotation and traction deformation, causing changes in the attitude (strike, dip, and dip angle) of the strata on both sides of the fault.

[0080] In application, the specific way to transform the construction topology rules into penalty functions is to introduce a step function or drift term related to the fault distance on both sides of the fault in the implicit equation for fault cutting relationships. A fault in three-dimensional space can be described by a surface equation, assuming this surface is Φ(x) = 0. Then, the position of any point in space relative to the fault can be determined by the sign of Φ(x). In traditional implicit modeling, the implicit function f(x) of the strata is continuous and smooth, but the existence of a fault breaks this continuity, and a step should be formed at the fault plane. By introducing a step function H(Φ(x)) or a drift term P(x) related to the fault distance, f(x) has different offset values ​​on both sides of the fault, expressing the displacement discontinuity of the strata on both sides of the fault.

[0081] In application, based on the movement characteristics of the fault, a trend rotation matrix is ​​added to the drift term to simulate the attitude rotation and traction deformation of the hanging wall relative to the footwall. In the dam site modeling example, the knowledge graph records the attitude and movement mode of the compressional-shear fault F1. Based on this information, the system adds a drift term containing a rotation matrix to the implicit equation. This matrix causes the hanging wall to rotate 15° relative to the footwall, a value derived from actual measurement data in geological surveys. When the implicit function is near the fault, the function value gradually changes along the fault strike direction according to the distance from the fault. Not only does it show a significant step at the fault plane, but the potential field values ​​of the strata on both sides of the fault also show a directional deflection, exhibiting the rotational deformation of the hanging wall strata relative to the footwall strata. The geological morphology around the fault modeled in this way can not only represent the basic characteristics of strata faulting but also reflect the secondary deformation phenomena accompanying fault activity.

[0082] Specifically, the complete process of transforming fault-cutting relationships into mathematical constraints is as follows: First, the system extracts the geometric features of the fault from the knowledge graph, including fault traces, fault attitude, and lithological combinations of the hanging wall and footwall. The system establishes a three-dimensional fault plane equation to determine the position and distance of any point in space relative to the fault. Second, the system determines the relative displacement direction of the hanging wall and footwall strata based on the direction and amount of fault movement. For example, for a normal fault, the hanging wall strata subside relative to the footwall strata, and the vertical position difference of the strata is the fault displacement. Third, the system adds fault constraints to the objective function of implicit modeling. When calculating the implicit function of the hanging wall rock mass, the implicit value near the fault is smaller than the implicit value far from the fault, forming a step. When calculating the implicit function of the footwall rock mass, the implicit value near the fault is larger than the implicit value far from the fault. This setting ensures that near the fault plane, the difference in implicit values ​​between the hanging wall and footwall strata is greater than 0, exhibiting a fault-related relationship. Fourth, the system adds a rotation matrix to the drift term in the implicit equation based on the movement characteristics of the fault. For compressional-torsional faults, the hanging wall strata will rotate clockwise or counterclockwise with the compressional-torsional movement. This rotation is encoded in the rotation matrix, so that the gradient direction of the implicit function on both sides of the fault will change accordingly. The resulting geological body morphology can reflect this rotational deformation.

[0083] The drift term can be expressed as a linear step function related to the fault distance, H(d) = fault displacement × sign(Φ(x)), where d is the distance from the point to the fault plane. Alternatively, a Gaussian smooth transition function can be used to ensure a smooth transition of the step near the fault plane rather than an abrupt jump. The rotation matrix can be determined by directly constructing it from the measured rotation angle parameters recorded in the knowledge graph, or by inverting the rotation angle from fault kinematics theory and borehole-measured attitude data to ensure the geological significance of the rotation parameters. The constraint penalty function can be implemented by calculating the difference between the implicit values ​​of the hanging wall and the implicit values ​​of the footwall at the sampling point; if the difference is negative, a secondary penalty is applied. Alternatively, a penalty can be applied to the continuity of the potential field values ​​on both sides of the fault, ensuring a smooth transition of the potential field gradient before and after the step, even though a step exists.

[0084] In terms of technical advantages, by explicitly encoding the geometric and kinematic features of faults in implicit equations, the model can accurately represent the cutting relationships of faults and associated tectonic deformations without erroneously generating coherent stratigraphic morphology on both sides of the fault. This is particularly important for complex tectonic regions with dense faults, preventing the model from artificially smoothing out fault-controlled tectonic features due to data fitting. Simultaneously, through measured rotational data recorded in the knowledge graph, the model can quantitatively reflect the associated deformations of fault activity, improving the prediction accuracy of stratigraphic morphology in tectonically complex areas and facilitating a more accurate assessment of fault risks in engineering projects.

[0085] As an optional implementation, the specific steps for assessing the confidence level of geophysical inversion data and converting it into soft constraints are as follows: In applications, geophysical inversion data refers to subsurface geological information obtained through physical exploration methods such as electromagnetic, seismic, and gravity surveys. Geophysical methods offer advantages such as large coverage area and high efficiency; however, their inversion results often exhibit ambiguity, meaning the same physical anomaly may correspond to multiple different geological models. Furthermore, geophysical data is susceptible to interference from factors such as overlying high-resistivity layers and noise, resulting in data quality less reliable than direct borehole sampling. Therefore, a confidence assessment is necessary before fusing geophysical data to identify high-confidence anomaly areas and eliminate low-confidence anomalies that may be noise. The key to this confidence assessment is utilizing the geophysical method applicability rules in the engineering geology knowledge graph. These rules define the applicable range of different geophysical methods under various geological conditions, as well as the mapping rules between physical properties and lithology, i.e., the range of physical properties that different lithologies and geological bodies should possess.

[0086] In application, the specific steps for confidence assessment of geophysical inversion data are as follows: First, query the geophysical method applicability rules and the mapping rules between physical properties and lithology in the engineering geology knowledge graph. The geophysical method applicability rules specify the applicable conditions for each geophysical method. For example, the CSAMT method is more applicable in areas with thick overburden, but less applicable in areas with strong magnetic interference. The physical property and lithology mapping rules specify the range of physical properties that different lithologies should possess. For example, the resistivity of granite is typically from several thousand Ω·m to tens of thousands of Ω·m. If CSAMT inversion shows that the resistivity of a certain area is only tens of Ω·m, then that area may correspond to low-resistivity lithologies such as fracture zones, mudstone, and karst, rather than granite. The system scores the confidence of each anomalous area in the geophysical inversion data based on two dimensions: first, the applicability of the geophysical method for that anomaly in that area; and second, the degree of matching between the lithological combination corresponding to the anomaly and the lithological combination that should exist in that stratum in the knowledge graph. In the dam site area embodiment, the system detected a low resistivity anomaly in the CSAMT (Constant Stress Mapping) system. The knowledge graph indicates that the bedrock in this area is dense limestone, which should have high resistivity, yet a low resistivity anomaly is observed. The system's query of the graph suggests that the low resistivity anomaly may correspond to a fault fracture zone or karst, and the knowledge graph confirms the presence of the F1 fault in this area. Combining this with the location of the fracture zone revealed by borehole ZK5, the system calculates that the spatial agreement between the low resistivity anomaly and the known fracture zone is 85%, and sets the confidence coefficient for this anomaly area at 0.8.

[0087] In application, the values ​​of retained high-confidence anomalies are normalized to the [0, 1] interval to characterize the probability of the corresponding geological attribute existing at that spatial location. For example, a geophysical anomaly area is identified as a fracture zone after confidence assessment, with a confidence coefficient of 0.8. The system directly uses this coefficient as the prior probability of the fracture zone existing at that spatial location, i.e., the probability of detecting a fracture zone at that location is 80%. Within the entire modeling domain, the system calculates the probability distribution of each spatial grid point being occupied by different geological attributes (such as intact rock mass, fracture zone, karst, etc.). This probability distribution is then incorporated into the Bayesian inference framework as prior information to guide the implicit modeling interpolation process.

[0088] Specifically, the complete workflow for geophysical data confidence assessment and probabilistic transformation is as follows: First, the system retrieves geological information about the modeling area from the knowledge graph, including lithology, occurrence, and known structures. Second, the system extracts features from each anomaly area retrieved through geophysical inversion, including the anomaly's location, shape, and intensity. Third, the system queries the applicability rules of the geophysical method to determine its effectiveness in the region. For example, CSAMT has poor applicability in areas with strong magnetic interference, and the system will lower the confidence level of the anomaly in that area. Fourth, the system queries the physical property and lithology mapping rules to determine whether the physical properties corresponding to the anomaly match the lithology that should exist in the region. For example, if the area should be hard rock but the inversion yields low resistivity, it is identified as a possible tectonic failure zone. Fifth, the system calculates the spatial correlation between the anomaly and known geological bodies (such as borehole points). If the spatial location of the anomaly and the specific geological body revealed by the borehole has a high degree of overlap, the confidence level is increased. In the dam site example, the CSAMT low resistivity anomaly matches the fractured zone revealed by borehole ZK5 with an 85% agreement, and the system sets the confidence coefficient to 0.8. Sixth, the system classifies the confidence levels of all anomaly regions, removing anomalies with confidence levels below a preset threshold (e.g., 0.5) and retaining high-confidence anomalies. Seventh, the system normalizes the confidence coefficients of the retained anomaly regions, mapping them to the [0, 1] interval to ensure their mathematical significance as a probability of existence as a geological attribute. Finally, the system inputs this probability distribution into the Bayesian inference framework, jointly driving the implicit modeling interpolation field update with borehole hard constraints.

[0089] The confidence score can be calculated by weighting the geophysical applicability score and the physical property / lithology matching score to obtain a comprehensive confidence coefficient. Alternatively, multiple scoring factors can be logically combined according to rules defined in the knowledge graph; for example, a high confidence score is obtained if both method applicability and physical property matching are satisfied, and a low confidence score is obtained if either is not satisfied. The anomaly threshold can be set using a unified global threshold based on the geological background and geophysical data noise level of the region, or by setting separate classification thresholds based on different lithologies and structural types, making the evaluation criteria more consistent with local geological characteristics. The probability distribution can be represented by directly converting high-confidence anomalies into the probability of the geological attributes present at that spatial location, or by establishing a Gaussian probability field, where the probability is highest at the anomaly center and gradually decreases outwards, forming a continuous probability distribution.

[0090] In terms of technical advantages, the ambiguity problem of geophysical data is systematically solved through knowledge graph-driven confidence assessment, avoiding arbitrary interpretations of geophysical anomalies. Simultaneously, by eliminating low-confidence anomalies, implicit modeling is not misled by geophysical noise, improving model reliability. By transforming high-confidence geophysical data into probability distributions and incorporating them into a Bayesian framework, the constraints of geophysical data are scientifically quantified. This effectively utilizes the coverage advantage of geophysical data, constrains data-free areas between boreholes, and prevents the excessive influence of low-quality geophysical data through the confidence mechanism.

[0091] As an optional embodiment, the implicit potential energy functional is configured as follows: In application, the implicit potential energy functional is the core optimization objective of implicit modeling, consisting of a weighted sum of three parts. The data fitting term constrains the implicit function to strictly pass through all borehole hard data points and approximates the probability distribution of the geophysical soft data in the sense of minimizing residuals. In the dam site example, there are 20 boreholes, each with clear lithology and stratigraphic information; these points are defined as hard constraints. The system ensures that the implicit function passes through every borehole, meaning the implicit function value is equal to the expected potential field value at that point. Simultaneously, after confidence assessment and probabilistic processing, the geophysical data becomes a spatially distributed probability field. The data fitting term approximates this probability field in the least-squares sense, ensuring that the spatial variation of the implicit function's gradient in the geophysical anomaly area matches the attribute variations reflected by the geophysical data.

[0092] In applications, the smoothing term, based on the magnitude of the implicit function's gradient or the Laplace operator, controls the geometric smoothness of the generated surface, preventing overfitting. The gradient of the implicit function represents the direction and rate of the fastest change in the implicit function value; a larger magnitude indicates a more drastic change in the implicit function at that point. Without a smoothing term, the optimization algorithm might generate a violently oscillating implicit function in an attempt to perfectly traverse all data points. This oscillation manifests as wrinkles and distortions on the interpolated geological interface, features with no geological significance and purely mathematical overfitting. The smoothing term constrains the smoothness of the generated surface by integrating the magnitude of the implicit function's gradient or integrating the Laplace operator (the sum of second-order partial derivatives). Weighting coefficients Control the intensity of smoothness. The larger the smoothing term, the stronger its influence and the smoother the generated surface. However, excessive smoothing can prevent the model from fitting the actual geological changes.

[0093] In application, the knowledge violation penalty term is calculated by deploying sampling points across the entire field or within a specific area of ​​interest, calculating the potential field deviation at each sampling point where the implicit potential field function violates stratigraphic sequence rules, morphological rules, or topological rules, and then summing all the violation deviations. This is the biggest innovation of this invention compared to traditional implicit modeling and is the core embodiment of knowledge-driven modeling. The system deploys a large number of sampling points uniformly or randomly within the modeling domain, potentially reaching millions. For each sampling point, the system checks whether the implicit function value at that point violates various geological rules. If it violates the sequence rule that Quaternary strata should overly overbedbed rock, the system calculates the deviation where the Q4 potential field value at that point is less than the T potential field value. This deviation is defined as max(0, T potential field value - Q4 potential field value). 2 The sum of the violation deviations of all sampling points is multiplied by the weighting coefficient. This results in a knowledge violation penalty. Weighting coefficient. The severity of penalties for violations of geological logic is controlled; the larger λ2 is, the stronger the penalty for violations.

[0094] Specifically, the complete mathematical expression of the energy functional is: in This is the data fitting term, mathematically equal to ‖Ax - b‖2 2 Here, matrix A is the design matrix composed of borehole coordinates and geophysical grid point coordinates, vector x is the coefficient vector of the implicit function, and vector b is the observation vector of the borehole and geophysical data. This term measures the overall fitting error between the implicit function and the observed data. For the smoothing term, mathematically it can be the integral of the magnitude of the implicit function gradient, i.e. It can also be in the form of the Laplace operator, i.e. . The knowledge violation penalty term is defined as the sum of the violation penalties over all sampled points, i.e. .

[0095] The error metric for the data fitting term can be the L2 form of Euclidean distance to measure the vertical deviation between the implicit function and the data points, or the Mahalanobis distance form to consider the anisotropic distribution of the data and assign different weights to data at different locations and directions. The smoothing term can use the magnitude of the gradient to constrain the rate of change of the implicit function, or the Laplace operator to constrain the curvature of the implicit function. The sampling strategy for knowledge violation penalties can be uniform sampling across the entire modeling domain to constrain the geological logic of the entire domain, or denser sampling in sparse borehole regions or regions of interest for more targeted constraints. The weight coefficients can be manually set based on data quality and the importance of geological rules, or automatically optimized using cross-validation or Bayesian methods to ensure a relatively balanced contribution from the three terms.

[0096] In terms of technical advantages, the energy functional consisting of the three components simultaneously optimizes three objectives: data fitting accuracy, surface smoothness, and geological logic compliance. Traditional implicit modeling only considers data fitting and smoothness, which can easily lead to geologically unreasonable patterns in sparse data areas. This method uses a knowledge violation penalty term to automatically guide the optimization process towards conforming to geological rules. By adjusting the relative magnitudes of the three weight coefficients, the balance between the three objectives can be flexibly controlled, and the modeling strategy can be adjusted according to different engineering application needs. For example, geological logic constraints can be strengthened for risk assessment, while data fitting can be strengthened for detailed exploration.

[0097] As an optional implementation, the iterative optimization algorithm employs a two-stage solution process using either the iterative reweighted least squares method or the penalty function method. In applications, iterative optimization algorithms for implicit potential energy functionals employ two mainstream methods: iterative reweighted least squares and the penalty function method. Both are standard solution algorithms for optimization problems containing inequality constraints or penalty terms. The basic idea of ​​iterative reweighted least squares is to transform constraints into weighted objective function terms, gradually strengthening the constraints by continuously adjusting the weights. The basic idea of ​​the penalty function method is to add penalty terms to the objective function, with the coefficient of the penalty term gradually increasing, making the violation of the original constraints progressively unacceptable.

[0098] In the application, the solution process is divided into two stages. In the initial iteration stage, the weight coefficient of the knowledge violation penalty term is set to a small value, allowing the algorithm to quickly obtain a geometrically approximate solution based on borehole and geophysical data. At this time, the primary goal of the optimization algorithm is to fit the hard data from the borehole and the soft data from the geophysical exploration, and quickly converge to obtain a 3D model with approximately correct geometry. In this stage, due to the knowledge violation penalty term... The initial solution is small, and the algorithm does not impose strong penalties for violations of geological logic rules, thus allowing for areas where the geological logic does not conform, such as dikes that may penetrate overburden. Although this initial solution is not perfect in terms of geological logic, it is basically correct in terms of geometric shape and data fitting accuracy, providing a good initial value for subsequent optimization.

[0099] In application, during subsequent iterations, the weight coefficient of the knowledge violation penalty term is gradually increased, forcing the implicit function to undergo asymptotic deformation while maintaining the data fitting quality, until the potential field bias of the semantic rule violation is reduced to below an acceptable threshold, and the total energy functional converges. The system may employ an exponentially increasing strategy to increase the weight coefficient. For example, multiplying by a factor of 1.5 in each iteration, or using a linear increment strategy, such as increasing by a fixed value in each iteration. As the energy increases, the penalty term energy for sampling points that violate geological rules increases rapidly. To minimize the total energy functional, the optimization algorithm is forced to adjust the parameters of the implicit function, causing changes in the implicit function values ​​at the violating sampling points and eliminating the violation. For example, for sampling points that violate sequence relations, the optimization algorithm will increase the potential field value of the overlying strata or decrease the potential field value of the underlying strata. This adjustment will gradually extend to nearby sampling points, forming an adjustment wave propagating from the local to the global level.

[0100] Specifically, taking the modeling of the dam site area as an example, the complete iterative optimization process is as follows. First, the system uses constraints from all borehole and geophysical data to set... =0, initiating the initial iteration. At this point, the energy functional does not contain a knowledge violation penalty term, and the optimization algorithm focuses on fitting the data. Due to the spatial inhomogeneity of the distribution of borehole and geophysical data, in some areas far from the data points, the value of the implicit function is determined by the shape of the interpolation kernel function. In this initial stage, the optimization algorithm iterates multiple times, and the energy gradually decreases, eventually converging to a local minimum. The generated model at this point is generally correct in shape, but there are geological logic violations. For example, in areas far from borehole ZK3, the implicit potential field value of the dike is larger than that of Q4, violating the rule that Q4 should overly overlay the dike. The system then begins the second stage of iteration based on this. The value is set to 0.01. The system starts again from the converged solution of the previous stage, using the same data constraints and smoothing terms, but with a smaller knowledge violation penalty term, for the next iteration. In this iteration, the system periodically evaluates the violation of geological rules at all sampling points. When a point is found where the potential field value of the dike is greater than the Q4 potential field value, the penalty term at that point generates an energy increment. To reduce the total energy, the optimization algorithm begins to adjust the potential field value of the dike at that point to decrease or the Q4 potential field value to increase. This adjustment gradually propagates outward, causing other sampling points near that point to gradually adjust their parameters, forming a local deformation wave centered on the violation area. After multiple iterations, the violation situation at all sampling points gradually improves. When the second stage is close to convergence, the system increases again. When the value reaches 0.05, the third iteration begins. At the start of each new stage, the optimization algorithm starts from the convergent solution of the previous stage and uses a stronger knowledge violation penalty to continue forcing the model to adjust towards geologically logical compliance. As the iteration progresses, As the number of violations gradually increases, the number of violations in the model gradually decreases, and the amount of violation bias gradually decreases. Ultimately, when... When the value is increased to a certain threshold (e.g., 100), the violation deviation at all sampling points is less than the preset acceptable threshold (e.g., 0.01). At this point, the model satisfies both the data fitting constraints and the geological logic constraints, the total energy functional converges to the minimum value, and the iterative optimization process is completed. In the final generated model, Q4 is clearly overlying the dike, which exhibits a plate-like morphology, the strata on both sides of the fault are correctly displaced, and all geological rules are satisfied.

[0101] The weighting strategy can employ exponential increments to rapidly increase the penalty intensity, promoting quick model deformation towards geologically compliant directions. Alternatively, a piecewise linear increment can be used to gradually increase the penalty intensity, avoiding excessive model deformation. The initial approximate solution can be obtained by using only borehole hard data to quickly obtain an interpolated solution centered on the borehole point. Alternatively, both borehole and geophysical data can be used simultaneously to obtain a broader-coverage initial solution. Convergence criteria can be based on the relative change of the total energy functional being below a preset threshold, such as 0.1%, or the maximum violation deviation of all sampling points being below a preset threshold, such as 0.01%. The iterative algorithm can be implemented using iterative reweighted least squares, with each iteration solving a weighted least squares problem. Alternatively, a penalty function method combined with the conjugate gradient method or other first-order optimization algorithms can be used to efficiently solve large-scale sparse linear equation systems.

[0102] In terms of technical advantages, the two-stage iterative strategy achieves a balance between data fitting and geological logic compliance. The initial stage quickly obtains a basically correct geometry, avoiding the undesirable local minima that might occur when directly solving a non-convex optimization problem with strong constraints. Subsequent stages gradually strengthen the constraints, making the model adjustment gradual and controllable, avoiding drastic non-physical oscillations. The advantages of this strategy lie in its good computational stability, fast convergence speed, and the fact that the final solution has both good data fitting accuracy and fully meets the geological logic requirements. Compared to solving large-scale optimization problems with strong constraints in a single step, this staged strategy has significant advantages in both computational complexity and convergence performance.

[0103] As an optional embodiment, the specific steps for constructing a hierarchical spatial index structure and performing local dynamic updates are as follows: In applications, constructing hierarchical spatial index structures is a key infrastructure for achieving efficient dynamic updates. Hierarchical spatial index structures, such as octrees, KD-trees, or multi-level voxel indexes, are data structures that recursively partition three-dimensional space. Octrees use a bipartitioning method to continuously divide the space into eight parts, with each node representing a cubic voxel and having eight child nodes representing the eight quadrants of that cube. KD-trees use an alternating binary partitioning method along different dimensions, with each node representing a hyperrectangle. Multi-level voxel indexes refer to multiple levels based on the size of the voxels, with the bottom layer being the finest voxel mesh and the upper layers being progressively coarser voxel aggregations. A common feature of these index structures is their ability to quickly locate the relationship between spatial locations and corresponding implicit equation nodes.

[0104] In the application, a fast retrieval mapping is established between spatial locations and corresponding implicit equation nodes. Implicit modeling typically employs the polynomial form of radial basis functions, where the implicit function can be represented as a weighted sum of a key node (called the center) and its corresponding coefficients. The computation of the implicit function value for the entire modeling domain depends on all these key nodes. Given a spatial location (x, y, z), calculating the implicit function value at that location requires traversing all key nodes, which is computationally intensive. Through hierarchical spatial indexing, it is possible to quickly determine which voxel a given spatial location falls into, thus requiring consideration only of the key nodes within that voxel and its neighboring voxels, significantly reducing computational load. In the dam site embodiment, the established octree index recursively divides the entire modeling domain into eight parts, starting from an initial maximum cube, until each minimum cube contains no more than a fixed number of key nodes (e.g., 10). This forms a coarse-to-fine tree structure, and the corresponding implicit equation nodes are organized according to this tree structure.

[0105] In application, when new exploration data is added to the database, the system uses an engineering geological knowledge graph to determine the type of logical conflict between the new data and the existing model. The system first reads the newly added borehole location and lithology information. Then, it compares the implicit function values ​​and corresponding geological bodies of the existing model at that borehole location with the lithology actually revealed by the new borehole. If the existing model predicts a dike but the new borehole reveals a cover layer, the system identifies this as a major logical conflict, potentially requiring a complete re-evaluation of the model. If the existing model predicts a dike with a thickness of 5m, but the new borehole reveals a dike with a thickness of 3m, the system identifies this as a geometric parameter correction conflict, which can be resolved through a local update.

[0106] In application, if the conflict is determined to be a geometric parameter correction type conflict, the affected spatial radius is determined with the new data coordinates as the center. The size of the radius depends on the influence range of the implicit function. For radial basis functions, the influence range is determined by the range of the kernel function and the anisotropy parameter. For example, if the effective influence range of the radial basis function at a certain point (defined as the kernel function value decaying to less than 5% of the peak value) is 30m, then the influence radius is set to 30m. In the dam site embodiment, the coordinates of the new borehole location are (K0+500, elevation 600). By calculating the influence range of the nearest critical node in the implicit function to this location, the influence radius is determined to be 30m.

[0107] In the application, an index structure is used to quickly retrieve all nodes falling within the influence radius. Starting from the root node of the octree, the system recursively determines whether the cube represented by each child node intersects with the sphere within the influence radius. If they intersect, the system recursively enters that child node and continues to check its child nodes. If they do not intersect, the system skips that child node and all its descendant nodes. This allows all relevant nodes to be found in logarithmic time complexity. In the dam site embodiment, octree retrieval quickly identifies all volume element nodes within a 30m radius around K0+500; the implicit equation coefficients corresponding to these nodes need to be recalculated.

[0108] In application, the implicit equations are reconstructed and solved only for the retrieved nodes, while keeping the potential field parameters of nodes outside the influence radius unchanged. The system first extracts all existing data points (mainly borehole points) and retrieved implicit equation nodes within the influence radius. Second, the system adds new data points to the constraints, forming a new set of constraint equations together with the existing data points. Third, the system only updates the implicit equation nodes within the influence radius, i.e., recalculating the coefficient values ​​corresponding to these nodes. Fourth, the system keeps the coefficients of all nodes outside the influence radius unchanged, meaning these nodes are completely unresponsive to the addition of new data. This achieves incremental local updates without needing to recalculate the implicit function of the entire field.

[0109] Specifically, the complete process of the local dynamic update is as follows. During construction, the TBM tunnel continuously passes through different geological layers during excavation. Each time it crosses an important geological interface, the surveyors record the lithology and depth at that location, forming new borehole data points. After the new data points are entered into the database, the system first checks the coordinates and lithology of the point and compares them with the existing 3D geological model. The system calculates the implicit function value of the existing model at the borehole location to infer the predicted lithology at that location. If the predicted lithology does not match the actual lithology, the system determines it as a conflict. The conflict type is analyzed; for example, a new borehole at K0+500 reveals an actual vein thickness of 3m, but the existing model predicts a thickness of 5m. This is a typical geometric parameter correction conflict. The system uses K0+500 as the center and, based on the range parameter of the vein's implicit function, determines the influence radius to be 30m. The system uses an octree index to quickly locate all voxel nodes within the radius. The system finds 200 key nodes within the influence range that need updating. The system extracts borehole data from these 200 nodes and their surrounding areas, adding newly added borehole points to form new constraint equations. The system then re-solves the implicit equations for these 200 nodes, calculating new coefficient values. The system keeps the coefficients of all nodes outside the influence radius unchanged. Finally, the system uses the updated implicit functions to re-extract the isosurfaces of the region, generating an updated 3D geological model. The entire update process takes 0.5 seconds, representing a leap from hours for a global model to seconds for a local update.

[0110] The index structure can be implemented using an octree, recursively traversing the tree nodes downwards to quickly determine the leaf nodes within the influence range. Alternatively, a KD-tree can be used, rapidly locating the spatial range by performing binary partitioning in different dimensions. The influence radius can be determined by setting a fixed radius based on the effective influence range of key nodes in the implicit function, or by adaptively calculating the actual influence range based on the location of new data, ensuring the radius size adapts to the curvature of the implicit function surrounding the data. The scope of resolving the equation can be determined by updating only the node coefficients within the influence radius, or by using locally weighted least squares to solve the equation for nodes within the influence radius, giving higher weights to nodes closer to the new data. Convergence can be determined by checking the change in the global energy functional; if the relative change is less than a preset threshold, convergence is considered achieved. Alternatively, the maximum change in the implicit function value at the updated node can be checked; if the maximum change is less than a preset threshold, convergence is considered achieved.

[0111] In terms of technological advantages, the update efficiency of 3D geological models has been fundamentally improved through a hierarchical spatial indexing structure and a local incremental reconstruction algorithm. Traditional methods must reconstruct and solve the entire implicit equation system after obtaining new data, resulting in a computational complexity of O(n log n). 3The computational complexity is significantly reduced by using an index structure for rapid localization, requiring only updates to local nodes within the influence radius. Simultaneously, it maintains unchanged parameters for nodes outside the influence radius, avoiding unnecessary computation of far-field nodes. This strategy enables even large-scale 3D geological models with millions of nodes to complete local corrections within seconds after acquiring new data, fully meeting the real-time update requirements during construction and providing technical support for the dynamic evolution of digital twin base plates.

[0112] Example 2 Based on the same inventive concept as Embodiment 1, Embodiment 2 of this application provides a geological implicit modeling and multi-source data dynamic interpolation system, which is used to implement the method described in Embodiment 1. Each functional module in Embodiment 2 corresponds one-to-one with the method steps in Embodiment 1: the data acquisition and knowledge graph management module carries the multi-source data acquisition and knowledge graph construction functions in Embodiment 1; the semantic parsing and mapping module implements the conversion function from hierarchical order, morphological, and structural rules to mathematical constraints in Embodiment 1; the multi-source data fusion module executes the Bayesian confidence assessment and probability field fusion functions in Embodiment 1; the energy functional construction and solution module undertakes the implicit potential field objective function construction and iterative optimization solution functions in Embodiment 1; and the local dynamic update module drives the hierarchical spatial index retrieval and local reconstruction within the ROI functions in Embodiment 1. This system solidifies the method steps of Embodiment 1 into a runnable computational architecture, with the specific module structure as follows: In application, a complete geological implicit modeling and multi-source data dynamic interpolation system constrained by a knowledge graph consists of five core functional modules working together. The data acquisition and knowledge graph management module is responsible for the system's input, collecting multi-source geological exploration data of the target modeling area. This data includes borehole data such as lithology, stratigraphy, and depth; geological profile data; and geophysical data such as CSAMT inversion results and seismic wave velocities. This module is also responsible for constructing or calling existing engineering geological knowledge graphs. Engineering geological knowledge graphs can be stored in graph databases such as Neo4j or in relational databases using RDF triples. The knowledge graph contains geological entities (strata, faults, dikes, etc.), geological attributes (age, occurrence, physical parameters, etc.), and semantic relationships between entities (overburden relationships, cutting relationships, contact relationships, etc.).

[0113] In the application, the semantic parsing and mapping module is responsible for parsing the geological semantic rules in the knowledge graph and transforming stratigraphic rules, geometric morphology rules, and tectonic topology rules into implicit potential field function boundary conditions, anisotropic parameters, and penalty functions, respectively. This module uses a rule engine or ontology inference engine to automatically query and reason about the knowledge graph. When the system receives a new modeling task, this module first retrieves all geological rules applicable to the region from the knowledge graph. For example, the query reveals that the strata in the region include three types: Q4, T, and γ, where Q4 overlies T and γ, and γ has a plate-like morphology. Next, the module processes these rules; for example, for overlying relation rules, it generates corresponding inequality constraints, and for morphological rules, it generates corresponding anisotropic matrix parameters.

[0114] In the application, the multi-source data fusion module is responsible for performing knowledge graph-based confidence assessment on the geophysical inversion data and fusing borehole hard constraints and geophysical soft constraints into a unified probabilistic interpolation field based on Bayesian inference theory. This module first calls the geophysical method applicability rules and property-lithology mapping rules generated by the semantic parsing and mapping module to score the confidence of each anomaly area in the geophysical inversion data. The confidence score considers multiple factors, including the applicability of the geophysical method in the location, the matching degree between the lithology corresponding to the anomaly and the expected lithology of the strata, and the spatial correlation between the anomaly and the borehole exposure results. This module removes anomalies with confidence scores below a threshold, retaining high-confidence anomalies. Then, the module normalizes the values ​​of high-confidence anomalies to the [0, 1] interval, representing the probability of the geological attribute's existence. Finally, this module uses the Bayesian inference framework to combine this probability distribution with the borehole hard constraints, outputting a fused probability field for use by the next module.

[0115] In the application, the energy functional construction and solution module is responsible for constructing an implicit potential field energy functional that includes data fitting terms, smoothing terms, and knowledge violation penalty terms. It then solves for the implicit function using an iterative optimization algorithm to extract isosurfaces and generate a 3D geological model. This module first receives data fitting constraints (probability fields of borehole hard data and geophysical soft data), smoothing parameters, and knowledge violation rules from the preceding modules. Based on these inputs, the module constructs the objective function: This module uses either iterative reweighted least squares or a penalty function method to solve the problem. The initial stage is set... The value is relatively small, and gradually increases in subsequent stages. The process continues until the total energy converges. After the solution is obtained, this module uses algorithms such as the moving cube to extract the isosurfaces of the implicit functions and generate a triangular mesh representation of the three-dimensional geological interface.

[0116] In the application, the local dynamic update module is responsible for maintaining a hierarchical spatial index structure covering the entire domain. When new exploration data is added, it quickly locates the affected areas and performs local incremental reconstruction, enabling real-time dynamic correction of the 3D geological model. After the initial model building is completed, this module constructs hierarchical index structures such as octrees, KD trees, or multi-level voxel indexes based on the key node distribution of implicit functions. This module is also equipped with a listener that is automatically triggered when new exploration data is added. The module first determines the conflict type between the new data and the existing model based on the logical rules of the knowledge graph. For conflicts involving geometric parameter correction, the module determines the radius of influence centered on the location of the new data, quickly retrieves all nodes within the radius using the index structure, resolves the implicit equations only for these nodes, keeps the external nodes unchanged, and completes the local incremental update.

[0117] Specifically, the collaborative operation of the five modules in handling a complete modeling project is as follows. After the dam site exploration project is launched, the data acquisition and knowledge graph management module first collects data from 20 boreholes and CSAMT geophysical data, and calls upon the existing engineering geological knowledge graph. The semantic parsing and mapping module is then activated, retrieving geological rules for the area from the knowledge graph, including Q4 overlying T and γ, γ exhibiting a plate-like morphology with NE30° SE dip angle of 75°, and F1 being a compression-torsional fault. This module transforms these rules into mathematical constraints; for example, the Q4 overlying rule is transformed into inequality constraints, the γ morphology rule into anisotropy matrices, and the fault rules into drift terms and rotation matrices. The multi-source data fusion module receives CSAMT low-resistivity anomaly data, performs confidence assessment using the geophysical rules of the knowledge graph, and finds that the low-resistivity anomaly has a high spatial agreement with the known fracture zone. The confidence coefficient is set to 0.8, and the anomaly is transformed into a probability field of fracture zone existence. The energy functional construction and solution module, integrating the above information, constructs an energy functional containing three terms and employs a two-stage iterative solution. The initial stage quickly obtains the model framework, while subsequent stages progressively strengthen geological constraints, ultimately generating a geologically logical 3D model. The local dynamic update module constructs an octree index after model completion. When the TBM reaches the K0+500 position and new dike thickness data is input, this module quickly locates the affected area and completes local model correction within 0.5 seconds. The collaboration of these five modules creates a complete closed loop from data input to model generation and dynamic correction, achieving automated, high-precision, and geologically logical 3D geological modeling and real-time updates.

[0118] The knowledge graph can be stored using graph databases such as Neo4j, which directly represent entities and relationships using nodes and edges, supporting complex graph queries. Alternatively, it can be stored in a relational database using RDF triples, with knowledge retrieval implemented through the SPARQL query language. The confidence assessment rules can be combined by weighted averaging of the scores of each assessment factor to obtain a comprehensive confidence score. Another approach is to judge multiple conditions in priority order, directly determining low confidence if any key condition is not met. The iterative optimization algorithm can employ iterative reweighted least squares, solving the problem using weighted least squares in each iteration. An augmented Lagrange method can also be used to dynamically balance constraints and the objective function. Isosurface extraction can utilize the moving cube algorithm to quickly generate triangular meshes, or it can use signed distance function rendering to directly generate smooth surfaces. The index structure can be maintained using static indexes, fixed after the initial model construction, or dynamic indexes, automatically adjusted after each local update.

[0119] In terms of technical advantages, the separation and collaboration of the five functional modules achieves a highly modular design, allowing each module to be developed, tested, and upgraded independently. Simultaneously, standardized interfaces between modules ensure excellent scalability, enabling flexible replacement or upgrades of individual modules without affecting others. For example, the implementation of the index structure can be easily replaced, or the performance of the knowledge reasoning engine can be upgraded. This design allows the system to adapt to different application needs and technological development trends. Furthermore, the modular design facilitates distributed parallel processing, allowing modules such as data fusion and energy optimization to be executed in parallel on multiple computers, improving overall system performance.

[0120] Example 3 Based on the same inventive concept as Embodiments 1 and 2, Embodiment 3 of this application provides a computer-readable storage medium on which a computer program is stored. When executed by a processor, the program implements the method described in Embodiment 1 or drives the system described in Embodiment 2 to run. Embodiment 3 encapsulates the method logic of Embodiment 1 and the system architecture of Embodiment 2 into an independent product form, facilitating direct deployment and application in digital platforms for various water conservancy and hydropower projects, transportation tunnels, and deep-earth mining projects. The specific implementation of this storage medium is as follows: In applications, computer-readable storage media refers to physical devices or storage media used to store computer programs, including various forms such as hard disk drives, solid-state drives, optical disks, magnetic tapes, and random access memory. When a computer program is executed by a processor, it transforms each step of the method described in any one of claims 1 to 8 into specific program instructions. These program instructions, through processor execution, invoke various hardware resources of the computer system, including the central processing unit (CPU) for logical operations and control, the graphics processing unit (GPU) for accelerating matrix calculations, memory for storing intermediate results, and the hard disk for persistently storing model data, thereby realizing the entire implicit geological modeling and dynamic interpolation method.

[0121] In application, the source code of this computer program can be written in programming languages ​​such as C, C++, Python, and Java. For the matrix operation-intensive optimization solution part, C++ or CUDA language is usually used to call the GPU for acceleration. For knowledge reasoning, data processing, etc., Python or Java is usually used, taking advantage of their rich third-party libraries. The overall structure of the program follows the five-module architecture of the system described in claim 9, and the computer program respectively implements functions such as data acquisition and knowledge graph management, semantic parsing and mapping, multi-source data fusion, energy functional construction and solution, and local dynamic update. During the execution of the program, firstly, the engineering geological knowledge graph data is loaded, and relevant rules are retrieved from the graph database through SPARQL or other query languages. Secondly, the program calls the semantic parsing module to generate implicit modeling constraints and parameters based on the retrieved rules. Thirdly, the program calls the multi-source data fusion module to evaluate the confidence of the geophysical data and fuse the borehole data. Then, the program constructs the energy functional and calls the optimization library (such as BFGS, LBFGS, or other nonlinear optimization algorithm libraries) to solve it. Finally, the program maintains the spatial index structure, listens for new data input, and performs partial updates.

[0122] Specifically, the execution flow of this computer program on a real computer system is as follows: Upon startup, initialization is performed first. The program reads the configuration file to obtain parameters such as the geographical extent and projected coordinate system of the modeling area. The program connects to the knowledge graph database and verifies its availability. The program allocates a memory buffer to store intermediate data such as borehole data, geophysical data, and coefficients of implicit functions. After initialization, the program enters the data loading stage. The program reads the coordinates, lithology, stratigraphy, and other information of 20 boreholes from the borehole data file and stores them in a data structure in memory. The program reads the resistivity value of each grid point from the geophysical data file (such as a CSAMT inversion grid file) and also stores it in memory. The program connects to the knowledge graph database through a database driver (such as graph-tool in Python or Neo4j driver in Java) and executes a SPARQL query to retrieve the geological rules of the area. Next, the program enters the semantic parsing stage. The program parses the retrieved stratigraphic relationship rules, such as...<Q4, overlies, T> The program generates corresponding inequality constraints. Based on morphological rules, it calculates the parameters of the anisotropic covariance matrix. Based on fault rules, it constructs fault drift terms and rotation matrices. Then, the program enters the data fusion stage. It calls a confidence assessment algorithm to score geophysical anomaly areas. Anomalies with confidence levels below a threshold are removed, and the remaining anomalies are normalized to generate a probability field. Based on Bayes' theorem, the program integrates borehole hard constraints and geophysical soft constraints. Next, the program enters the energy functional construction and solution stage. Based on the preceding constraints and parameters, the program constructs the objective function E(f). The program initializes the weight coefficients. and The program initializes the objective function. It then calls an optimization library (such as the `minimize` function in SciPy or a self-implemented iterative reweighted least squares algorithm) to perform iterative solutions. In each iteration, the program calculates the gradient of the objective function and updates the coefficients of the implicit function. The program periodically evaluates violations of geological rules at sampling points and calculates the value of the knowledge violation penalty term. Based on the set convergence criteria, the program determines whether to increase the value of the penalty term. Moving to the next stage. After convergence, the program uses the moving cube algorithm to extract isosurfaces, generating a 3D geological interface represented by a triangular mesh. Finally, the program enters the dynamic update stage. The program constructs an octree index based on the key node distribution of implicit functions. The program starts a listening thread to monitor for new data input. When new borehole data arrives, the program determines the conflict type, identifies the radius of influence, quickly retrieves relevant nodes using the index, and resolves only for these nodes, completing the local update. The program outputs the updated model to a file or database for subsequent engineering applications.

[0123] The program can be executed sequentially using a single thread, following the order of data loading, semantic parsing, data fusion, energy optimization, and dynamic updates. Alternatively, it can be executed in parallel using multiple threads, allowing independent modules such as data loading, knowledge graph querying, and energy optimization to run concurrently on different processor cores, improving overall execution efficiency. Intermediate data can be stored in memory, saving all intermediate results in RAM to accelerate the computation process. A hybrid memory-disk storage approach can also be used, storing some intermediate results on disk for very large-scale models while keeping memory usage within a reasonable range. The program's output format can use standard 3D model formats such as STL or OBJ for easy viewing in common 3D visualization tools, or a proprietary binary format for efficient storage and fast loading. Optimization methods can utilize multi-core parallel optimization libraries on the CPU to fully leverage the computing power of modern multi-core processors, or GPU acceleration using CUDA or OpenCL for large-scale matrix operations and optimization calculations on the graphics processing unit.

[0124] In terms of technological advantages, this computer program fully implements complex geological modeling methods into an executable software product, enabling the method to run on general-purpose computer hardware. The program's modular design allows for independent debugging and optimization of each functional component. It supports combinations of multiple programming languages, leveraging the efficiency of C++ for numerical computation and the flexibility of Python for data processing and knowledge reasoning. The program's parallel processing capabilities allow even large-scale modeling problems containing millions of implicit function key nodes to be completed within a reasonable timeframe on modern multi-core or many-core computers. The program's dynamic update mechanism allows for the updating of the entire 3D geological model within seconds of the addition of new data during construction, meeting the needs of real-time engineering decision-making. Furthermore, the program's readable storage medium format allows for the widespread dissemination and application of this method, unrestricted by specific software platforms, which is beneficial for the industrialization and intellectual property protection of this technology.

[0125] The above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application, and should all be included within the protection scope of this application.

Claims

1. A method for implicit geological modeling and dynamic interpolation of multi-source data, characterized in that, Includes the following steps: Acquire multi-source geological exploration data of the target area, and construct or call an engineering geological knowledge graph, which covers geological entities, geological attributes, and semantic relationships between entities; The geological semantic rules in the engineering geology knowledge graph are analyzed, and the mapping relationship between semantic rules and mathematical constraints is established. The sequence relationship rules, geometric morphology rules and structural topology rules are transformed into the boundary conditions, anisotropic parameters and penalty functions of the implicit potential field function, respectively. The multi-source geological exploration data is fused based on Bayesian inference theory. The borehole exploration data is defined as deterministic hard constraints, and the geophysical inversion data is transformed into probabilistic soft constraints after confidence evaluation using the applicability rules in the engineering geology knowledge graph. An implicit potential field energy functional is constructed, which includes a data fitting term, a smoothing term, and a knowledge violation penalty term. The knowledge violation penalty term is used to measure the degree to which the implicit potential field function violates the geological semantic rules. The energy functional is solved by an iterative optimization algorithm to obtain the implicit function expression of the geological interface. Based on the solved implicit function, isosurfaces are extracted to generate a three-dimensional geological model, and a hierarchical spatial index structure is constructed. The three-dimensional geological model is then dynamically updated locally based on newly added exploration data.

2. The method according to claim 1, characterized in that, The specific method for transforming hierarchical relationship rules into boundary conditions is as follows: If the engineering geological knowledge graph defines stratum A as overlying stratum B, then an inequality constraint is constructed in the normal direction of the stratum contact surface, which forces the implicit potential field value of stratum A at any spatial point to be greater than the potential field value of stratum B, in order to prevent non-geologically induced stratum inversion phenomena from occurring during the modeling process.

3. The method according to claim 1, characterized in that, The specific method for converting geometric shape rules into anisotropic parameters is as follows: Based on the geological body genetic types and occurrence attributes recorded in the engineering geological knowledge graph, including strike, dip and dip angle, an anisotropic covariance matrix or weight matrix is ​​constructed. The matrix is ​​used to control the range ratio and rotation angle of the interpolation kernel function in different spatial directions so that the generated geological body morphology conforms to its sedimentary or tectonic background. Specifically, for sedimentary geological bodies extending in a layered manner, the matrix makes the range of the interpolation kernel function along the strike-dip direction of the bedding plane significantly greater than the range along the normal direction; for intrusive geological bodies occurring in vein-like patterns, the matrix makes the interpolation kernel function flattened into an ellipsoid, consistent with the occurrence of the vein.

4. The method according to claim 1, characterized in that, The specific method for transforming the topology rules into penalty functions is as follows: For fault cutting relationships, a step function or drift term related to the fault distance is introduced into the implicit equation on both sides of the fault to express the displacement discontinuity of the strata on both sides of the fault. Based on the motion properties of the fault, a trend rotation matrix is ​​added to the drift term to simulate the attitude rotation and traction deformation of the hanging wall relative to the footwall, thereby achieving a complete mathematical expression of the fault cutting relationship and associated structural deformation.

5. The method according to claim 1, characterized in that, The specific steps for assessing the confidence level of geophysical inversion data and converting it into soft constraints are as follows: Query the geophysical method applicability rules and the mapping rules between physical properties and lithology in the engineering geology knowledge graph, and score the confidence level of each anomaly area in the geophysical inversion data. Remove noise anomalies with a confidence level below a preset threshold; The values ​​of the retained high-confidence anomaly areas are normalized to the range of zero to one to characterize the probability of the corresponding geological attributes existing in that spatial location; The normalized probability distribution mentioned above is introduced as prior information into the Bayesian inference framework and, together with the borehole hard constraint, drives the interpolation field update of the implicit model, thereby using geophysical data to constrain model inference for areas without borehole coverage.

6. The method according to claim 2, 3 or 4, characterized in that, The implicit potential energy functional is composed of a weighted sum of three parts: The data fitting term constrains the implicit function to strictly pass through all borehole hard data points and approximates the probability distribution of geophysical soft data in the sense of minimizing residuals; The smoothing term, based on the modulus of the implicit function gradient or the Laplacian operator, controls the geometric smoothness of the generated surface and prevents overfitting. The knowledge violation penalty item is obtained by setting up sampling points in the entire field or a specific area of ​​interest, calculating the potential field deviation of the implicit potential field function at each sampling point that violates the hierarchical rule, the morphological rule, or the topological rule, and then summing all the violation deviations by weight. The energy functional achieves a coordinated balance between data fitting accuracy, surface smoothness, and geological logic compliance by adjusting the weight coefficients of the smoothing term and the knowledge violation penalty term.

7. The method according to claim 6, characterized in that, The iterative optimization algorithm employs either iterative reweighted least squares or a penalty function method, and the solution process is divided into two stages: In the initial iteration phase, the weight coefficient of the knowledge violation penalty item is set to a small value, so that the algorithm can quickly obtain a geometric approximate solution based on borehole and geophysical data. In subsequent iterations, the weight coefficient of the knowledge violation penalty term is gradually increased, forcing the implicit function to undergo gradual deformation while maintaining the data fitting quality, until the potential field deviation that violates the semantic rules is reduced to below an acceptable threshold, the total energy functional converges, and finally the implicit function solution that satisfies both data constraints and geological logic constraints is output.

8. The method according to claim 1, characterized in that, The specific steps for constructing the hierarchical spatial index structure and performing local dynamic updates are as follows: Pre-construct an octree, KD tree, or multi-level volume element index structure for the entire modeling domain to establish a fast retrieval mapping between spatial locations and corresponding implicit equation nodes; When new exploration data is entered into the database, the type of logical conflict between the new data and the existing model is determined based on the engineering geological knowledge graph. If the conflict is determined to be a geometric parameter correction type conflict, the affected spatial radius is determined with the new data coordinates as the center, and the index structure is used to quickly retrieve all nodes that fall within the affected radius. The implicit equations are reconstructed and solved only for the nodes retrieved above, while keeping the potential field parameters of the nodes outside the influence radius unchanged; Through the aforementioned incremental local reconstruction, the three-dimensional geological model can be rapidly corrected without requiring a full recalculation.

9. A geological implicit modeling and multi-source data dynamic interpolation system under knowledge graph constraints, characterized in that, include: The data acquisition and knowledge graph management module is used to collect multi-source geological exploration data of the target area and construct or call up engineering geological knowledge graphs; The semantic parsing and mapping module is used to parse the geological semantic rules in the engineering geological knowledge graph, and to transform the sequence relationship rules, geometric morphology rules and structural topology rules into the boundary conditions, anisotropy parameters and penalty functions of the implicit potential field function, respectively. The multi-source data fusion module is used to perform knowledge graph-based confidence assessment on geophysical inversion data and to fuse borehole hard constraints and geophysical soft constraints into a unified probabilistic interpolation field based on Bayesian inference theory. The energy functional construction and solution module is used to construct an implicit potential energy functional that includes data fitting terms, smoothing terms, and knowledge violation penalty terms, and solve the implicit function through an iterative optimization algorithm to extract isosurfaces and generate a three-dimensional geological model. The local dynamic update module is used to maintain a hierarchical spatial index structure covering the entire area. When new exploration data is added to the database, the affected area is quickly located and local incremental reconstruction is performed to realize the real-time dynamic correction of the three-dimensional geological model. The data acquisition and knowledge graph management module, semantic parsing and mapping module, multi-source data fusion module, energy functional construction and solution module, and local dynamic update module work together to execute the method of any one of claims 1 to 8.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method of any one of claims 1 to 8, or drives the system of claim 9 to operate.

Images