Deep phosphor ore tectonic stratum three-dimensional prediction system based on lithofacies paleogeographic constraint

Abstract

The application discloses a deep phosphorite tectonic stratum three-dimensional prediction system based on lithofacies paleogeographic constraint, belongs to the three-dimensional prediction technical field, and comprises a data input module, a stratum analysis module, a prediction modeling module and a delineation module. With the depth of a sequence interface as the center and after Gauss attenuation is performed on a Gauss attenuation width determined by the thickness difference of upper and lower strata, the coded values of the attenuated sedimentary facies types are weighted and superimposed with the contribution weight of each sequence interface, and then the thickness of the stratum is taken as an auxiliary constraint term for weighted summation to establish a three-dimensional space paleogeographic constraint field. The prediction modeling module constructs a three-dimensional prediction model in three nonlinear modes of power index scaling, fracture index modulation and direct coupling of the three-dimensional space paleogeographic constraint field. The delineation module performs boundary delineation through an anisotropic three-dimensional variation function model to generate a closed three-dimensional ore body envelope surface, solves the problem that the paleogeographic constraint is disconnected with the three-dimensional prediction, and realizes fine three-dimensional prediction of deep phosphorite tectonic strata.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional prediction technology, and in particular to a three-dimensional prediction system for deep phosphate rock tectonic strata based on lithofacies paleogeographic constraints. Background Technology

[0002] Phosphate rock is an important strategic mineral resource in my country, widely used in phosphate fertilizer production and fine phosphate chemical industry. After years of large-scale mining, shallow phosphate rock resources are becoming increasingly depleted, and deep, concealed phosphate rock has become the main focus of current exploration work. Deep phosphate rock has complex structures and varied spatial morphology, and traditional prediction methods based on two-dimensional profiles (such as the two-section method and the bank slope structure method) are difficult to meet the spatial positioning accuracy requirements of deep mineral exploration. Therefore, three-dimensional prediction methods are gradually becoming the research focus in this field.

[0003] Existing three-dimensional prediction methods for deep phosphate deposits typically include three stages: three-dimensional geological model construction, extraction of ore-controlling elements, and mineralization prediction. Regarding three-dimensional geological model construction, existing technologies utilize borehole data and geophysical data to establish a three-dimensional geological model, and then extract geological information through spatial interpolation or geostatistical methods. In terms of mineralization prediction, existing technologies combine machine learning algorithms with mineral exploration theories to construct a comprehensive three-dimensional prediction model.

[0004] However, the formation of sedimentary phosphate deposits is strictly controlled by the paleogeographic environment, and the distribution, thickness, and grade of ore-bearing rock strata are closely related to specific sedimentary facies zones. Sequence stratigraphy and lithofacies paleogeographic patterns are key geological factors controlling the spatial location of phosphate deposits, but in existing three-dimensional prediction methods, paleogeographic analysis typically remains at the qualitative descriptive level, and the geological constraints it contains cannot be used in three-dimensional spatial calculations. Specifically: First, the stratigraphic unit division in three-dimensional geological models mainly relies on borehole interpolation, and the spatial location of sequence boundaries is not constrained and corrected by the paleogeographic pattern, resulting in insufficient accuracy of the stratigraphic model; Second, subsequent three-dimensional mineralization prediction models are disconnected from the paleogeographic analysis results, with paleogeographic information only used for qualitative judgment of exploration direction in the early stages and not participating in the calculation of prediction index values. Summary of the Invention

[0005] In response to the above situation, this invention establishes a paleogeographic constraint field that can participate in three-dimensional spatial calculations, enabling paleogeographic constraint information to permeate the entire process of stratigraphic division and mineralization prediction, thereby achieving the integration of paleogeographic constraints and three-dimensional prediction and improving the accuracy and reliability of three-dimensional prediction of deep phosphate ore tectonic strata.

[0006] The technical solution is a three-dimensional prediction system for deep phosphate rock tectonic strata based on lithofacies paleogeographic constraints, including: The data input module acquires multi-source geological data for the target area, including borehole lithology data, well logging curve data, seismic reflection data, gravity data, and magnetic data. The stratigraphic analysis module identifies the sequence boundaries of the target area based on the multi-source geological data. It performs Gaussian attenuation with the depth of each sequence boundary as the center and the Gaussian attenuation width determined by the difference in strata thickness above and below each sequence boundary. The attenuated sedimentary facies type codes are then weighted and superimposed with the contribution weight of each sequence boundary, and finally weighted and summed with strata thickness as an auxiliary constraint term to establish a three-dimensional paleogeographic constraint field. The stratigraphic analysis module also transforms stratigraphic boundary identification into solving the local maximum surface of the three-dimensional paleogeographic constraint field along the depth direction. That is, the depth points at each horizontal position in the three-dimensional space where the three-dimensional paleogeographic constraint field reaches a local maximum along the depth direction are used to form the stratigraphic unit boundary surface, and the top and bottom depths of each stratigraphic unit are determined according to the stratigraphic unit boundary surface. The prediction modeling module fits model parameters based on known phosphate deposit control parameters and grade observation values, extracts phosphate control parameters for deep phosphate deposits in the target area, and generates three-dimensional distribution prediction results for phosphate resources through a nonlinear prediction model that includes power scaling, fracture index modulation, and direct coupling of three-dimensional paleogeographic constraint fields. The delineation module establishes an anisotropic three-dimensional variogram model and calculates the predicted value of the target location using Kriging interpolation based on the anisotropic three-dimensional variogram model. The predicted value of the target location is compared with the industrial grade index value of phosphate rock. Spatial locations that are greater than or equal to the industrial grade index value of phosphate rock are marked as the predicted phosphate ore body boundary, and a closed three-dimensional ore body envelope surface is generated.

[0007] Furthermore, the expression for the three-dimensional spatial paleogeographic constraint field is: ; in, , , These represent the horizontal x-coordinate, horizontal y-coordinate, and continuous depth coordinates in three-dimensional space, respectively. This refers to the number of the layer sequence interface; This represents the total number of hierarchical interfaces. For three-dimensional paleogeographic constraint fields in coordinates The value at; For the first The depth constant value of each hierarchical interface; For the first The Gaussian decay width of each sequence interface, and . ,in For the first The thickness difference between the strata above and below the sequence boundary. For all The largest stratigraphic thickness difference at each sequence boundary. The baseline Gaussian attenuation width, For adaptive adjustment index and ; For the first The contribution weights of each sequence interface; exp(·) is a function of natural constants. An exponential function with base 0; coordinates ( The sedimentary facies type code value at () is shown. coordinates The thickness of the strata at that location; This is an auxiliary weighting coefficient for formation thickness; The continuous depth coordinates With respect to the constant value of the layered interface depth For depth parameters with the same dimensions, when the continuous depth coordinates The value is equal to the constant value of the interface depth of a certain sequence. At that time, the Gaussian decay term corresponding to the sequence interface in the three-dimensional paleogeographic constraint field reaches its maximum value.

[0008] Furthermore, the formation analysis module includes: The sequence interface identification unit calculates the difference between the natural gamma logging value and the mean of the sliding window, and the difference between the density logging value and the mean of the sliding window at each depth point based on the well logging data and the seismic reflection data. The two values ​​are multiplied and the weighted average is taken within the sliding window as the sequence interface identification response value. The local peak position of the sequence interface identification response value is taken as the sequence interface position. The sedimentary facies analysis unit performs sedimentary facies analysis on each sequence unit based on the borehole lithology data and the identified sequence boundaries, and generates sedimentary facies type coding values. The paleogeographic pattern restoration unit, based on the sequence boundary identification results, sedimentary facies type coding values, and stratigraphic thickness, performs Gaussian attenuation with the depth of each sequence boundary as the center and the Gaussian attenuation width determined by the difference in stratigraphic thickness above and below each sequence boundary. The attenuated sedimentary facies type coding values ​​are then weighted and superimposed with the contribution weight of each sequence boundary, and finally weighted and summed with stratigraphic thickness as an auxiliary constraint term to establish the three-dimensional spatial paleogeographic constraint field.

[0009] Furthermore, for sedimentary phosphate deposits in the target area, the mapping relationship of the sedimentary facies type coding values ​​is as follows: intertidal facies is coded as 5, subtidal facies as 4, platform facies as 3, slope facies as 2, and basin facies as 1.

[0010] Furthermore, the expression for the three-dimensional distribution prediction result of the phosphate rock resources is as follows: ; in, , , These represent the horizontal x-coordinate, horizontal y-coordinate, and depth coordinates in three-dimensional space, respectively. coordinates Three-dimensional distribution prediction of phosphate rock resources at the location; Numbering of ore-controlling parameters; The total number of ore-controlling parameters; For the first Individual mineral control parameters in coordinates The normalized value at that point; coordinates Fault structure density at the location; For the first The power-law adjustment coefficients of the mineral control parameters; exp(·) is an exponential function with the natural constant e as the base; For the first The fracture index adjustment coefficient of individual mineral control parameters; For the first The fusion weight of individual mineral control element parameters; The coupling coefficient of the three-dimensional paleogeographic constraint field is denoted as . The three-dimensional paleogeographic constraint field in coordinates The value at the location; the coupling coefficient satisfy ,in For the first Correlation coefficients between individual mineral control parameters and the three-dimensional paleogeographic constraint field; The power index adjustment coefficient The range of values ​​is and The fracture index adjustment coefficient The value of is a non-negative number; the fusion weight satisfy And each All are greater than or equal to 0.

[0011] Furthermore, the ore-controlling parameters extracted by the prediction modeling module include the thickness of the ore-bearing rock series, the thickness of the phosphate rock layer, the grade of the phosphate rock layer, the dip angle of the strata, the density of fault structures, and the distance of paleogeographic facies zones. The fracture structure density parameter is determined by the three-dimensional Gaussian kernel fracture density calculation, that is, each fracture is taken as a three-dimensional Gaussian distribution source, a three-dimensional Gaussian kernel is established with the fracture spatial coordinates as the center, and all Gaussian kernels are superimposed and divided by the Gaussian kernel volume to obtain the three-dimensional continuous distribution of fracture structure density. The paleogeographic facies distance parameter is determined by weighted facies distance calculation, that is, for each sedimentary facies type, the shortest planar distance from the target location to the boundary of that sedimentary facies type is calculated, and the distance is summed by weighting the ore-controlling weight coefficients of each facies type.

[0012] Furthermore, the prediction modeling module also includes a model iteration unit, which iteratively optimizes the joint optimization objective function using the alternating direction multiplier method. In each iteration, the fusion weight, the power exponent adjustment coefficient, and the fracture exponent adjustment coefficient are updated sequentially. The fusion weight update uses the proximal gradient operator to ensure its non-negativity and the constraint that the summation is 1.

[0013] 8. The system according to claim 2, wherein the expression of the anisotropic three-dimensional variogram model is: ; in, It is a three-dimensional spatial lag distance vector. , , These represent the components of the lag distance vector in three-dimensional space along the horizontal x-axis, horizontal y-axis, and depth directions, respectively. It belongs to the same depth dimension as the depth coordinate Z; The hysteresis distance vector The corresponding variogram value; This is the base value; For variable range; It is an anisotropic tensor matrix; For the natural constant An exponential function with base 0; The hysteresis distance vector Transpose of; It is a quadratic form of the hysteresis distance vector and the anisotropic tensor matrix; The anisotropic tensor matrix The orthogonal rotation matrix is ​​obtained by eigenvalue decomposition of the diagonal eigenvalue matrix and the orthogonal rotation matrix, wherein the orthogonal rotation matrix is ​​the gradient vector field formed by taking the partial derivatives of the three-dimensional paleogeographic constraint field along each coordinate axis. Principal component analysis is used to determine the following: the three components of the gradient vector are used as three feature dimensions. A data matrix is ​​constructed from the three component values ​​of all grid points, and principal component analysis is performed to obtain the three principal component directions as the three column vectors of the orthogonal rotation matrix. The three diagonal elements of the diagonal eigenvalue matrix are the reciprocals of the squares of the variable range parameters of the three principal directions. in = ( ).

[0014] Furthermore, it also includes: The 3D display module generates a 3D geological model and a 3D phosphate resource distribution prediction map of the target area based on the 3D ore body envelope. The prediction and verification module determines the preferred exploration target area based on the predicted three-dimensional distribution of the phosphate rock resources and the three-dimensional ore body envelope, and outputs the three-dimensional spatial coordinate range of the preferred exploration target area.

[0015] Furthermore, it includes the following steps: Step 1: Obtain multi-source geological data for the target area, including borehole lithology data, well logging curve data, seismic reflection data, gravity data, and magnetic data; Step 2: Identify sequence boundaries based on the multi-source geological data. Using the depth of each sequence boundary as the center and the Gaussian attenuation width determined by the difference in strata thickness above and below each sequence boundary, perform Gaussian attenuation. Then, weight the attenuated sedimentary facies type codes with the contribution weight of each sequence boundary and sum them up with strata thickness as an auxiliary constraint term to establish a three-dimensional paleogeographic constraint field. The local maxima of the three-dimensional paleogeographic constraint field along the depth direction constitute the stratigraphic unit boundary surface, and the top and bottom depths of each stratigraphic unit are delineated accordingly. Step 3: Fit model parameters based on known ore-controlling parameters and grade observation values ​​of phosphate deposits. Use the normalized values ​​of each ore-controlling parameter as the base and the power-law adjustment coefficient as the exponent for power-law scaling. Use the exponential function value of the fracture structure density for fracture-law modulation. Then, directly superimpose the three-dimensional paleogeographic constraint field into the prediction results to establish a three-dimensional prediction model and generate a three-dimensional distribution prediction result of phosphate resources. Step 4: Establish an anisotropic three-dimensional variogram model, calculate the predicted value of the target location through Kriging interpolation, compare the predicted value of the target location with the industrial grade index value of phosphate rock, mark the spatial location that is greater than or equal to the industrial grade index value of phosphate rock as the predicted boundary of the phosphate ore body, and generate a closed three-dimensional ore body envelope surface.

[0016] Due to the adoption of the above technical solutions, the present invention has the following advantages compared with the prior art; 1. Using the depth of each sequence boundary as the center and the Gaussian attenuation width determined by the difference in strata thickness above and below the sequence boundary as the Gaussian attenuation width, the attenuated sedimentary facies type codes are weighted and superimposed with the contribution weight of each sequence boundary, and then weighted and summed with strata thickness as an auxiliary constraint term to establish a three-dimensional paleogeographic constraint field. The magnitude of this constraint field directly reflects the comprehensive control intensity of paleogeographic conditions on deep strata. At the same time, the stratigraphic boundary identification is transformed into solving the local maximum surface of this constraint field along the depth direction, so that paleogeographic constraints directly drive stratigraphic division. In the prediction stage, the paleogeographic constraint field is directly superimposed into the prediction results as an independent prediction term, realizing the direct participation of paleogeographic constraint information in the entire process from stratigraphic division to mineralization prediction, and upgrading paleogeographic analysis from qualitative description to computable numerical constraints.

[0017] 2. Power-law scaling applies adaptive nonlinear transformations to ore-controlling elements using learnable power-law exponents. When the power-law exponent is greater than 1, it amplifies high-value regions; when it is less than 1, it compresses low-value regions, enabling the model to characterize the threshold response characteristics of different ore-controlling elements. Fault exponent modulation modulates the contribution of each ore-controlling element with the exponential function value of fault density, making the promoting effect of fault development on mineralization grow exponentially, which matches the principle of mineralization kinetics. Direct coupling of three-dimensional paleogeographic constraint fields directly superimposes the paleogeographic constraint fields in their original numerical form, avoiding the weakening of their physical meaning in nonlinear transformations. This enables the prediction model to adaptively characterize complex nonlinear mineralization laws. Attached Figure Description

[0018] Figure 1 This is a flowchart of the system modules of the present invention. Figure 2 This is a flowchart of the formation analysis module of the present invention.

[0019] Figure 3 This is a diagram illustrating the method steps of the present invention. Detailed Implementation

[0020] The foregoing and other technical contents, features and effects of the present invention are described in conjunction with the appendix below. Figures 1 to 3 The detailed description of the embodiments will make this clear. All structural details mentioned in the following embodiments are based on the accompanying drawings.

[0021] Based on existing technology, in Example 1, this example provides a three-dimensional prediction system for deep phosphate rock tectonic strata based on lithofacies paleogeographic constraints, including a data input module, a stratigraphic analysis module, a prediction modeling module, and a delineation module.

[0022] The data input module is used to acquire multi-source geological data for the target area. This multi-source data includes borehole lithology data, well logging data, seismic reflection data, gravity data, and magnetic data. Borehole lithology data includes lithology name, lithology code, core recovery rate, and rock color description, used to identify stratigraphic lithological assemblages and sedimentary facies types. Lithology name and lithology code directly reflect the characteristics of rock assemblages in different sedimentary environments, while core recovery rate is used to assess the representativeness and reliability of borehole data. Well logging data includes natural gamma ray logging, resistivity logging, sonic transit time logging, and density logging. Natural gamma ray logging reflects changes in clay content, while density logging reflects changes in rock density. Both resistivity and sonic transit time logging (STL) data are sensitive to sequence interfaces. Resistivity logging curves and sonic transit time logging curves are used to help identify lithological interfaces and changes in fluid properties. Seismic reflection data includes seismic reflection time profiles and seismic wave velocity data. Seismic reflection time profiles provide the continuous reflection structure of subsurface strata, while seismic wave velocity data is used for time-depth conversion and lithology prediction. Gravity data includes Bouguer gravity anomaly data and free-air gravity anomaly data, which are used to identify regional fault structures and basement undulations. Magnetic data includes magnetic anomaly data and magnetic susceptibility data, which are used to identify the spatial distribution of concealed rock masses and fault structures. The above multi-source geological data provide the data foundation for three-dimensional spatial geological modeling and paleogeographic constraint field construction in this system.

[0023] The stratigraphic analysis module is used to identify sequence boundaries in the target area based on the multi-source geological data. It performs Gaussian attenuation centered on the depth of each sequence boundary and using a Gaussian attenuation width determined by the difference in strata thickness above and below each sequence boundary. The attenuated sedimentary facies type codes are then weighted and summed according to the contribution weight of each sequence boundary, and finally weighted and summed again using strata thickness as an auxiliary constraint term to establish a three-dimensional paleogeographic constraint field. The physical meaning of this three-dimensional paleogeographic constraint field is the comprehensive control intensity of the sedimentary environment on deep strata. Its numerical value reflects the degree to which a specific spatial location is controlled by paleogeographic conditions. High-value areas correspond to favorable sedimentary mineralization environments, while low-value areas correspond to unfavorable environments. This elevates paleogeographic analysis from qualitative description to a numerical field that can participate in three-dimensional spatial calculations, solving the problem that paleogeographic constraint information cannot be integrated into quantitative prediction models in traditional methods. In sequence stratigraphy, sequence boundaries are interfaces where sedimentary environments undergo abrupt changes. Their control intensity over underlying strata decreases with increasing depth. This system uses a Gaussian function to characterize this attenuation process because the Gaussian function... Gaussian decay exhibits characteristics of smooth decay, infinite differentiability, and clear physical meaning of parameters. Its decay pattern matches the geological law that the control intensity of sedimentary interfaces smoothly decreases with increasing distance. Compared with linear decay or step functions, Gaussian decay better reflects the continuity and gradualness of sedimentation. The contribution weight of each sequence boundary is determined by normalizing the significance score of each sequence boundary using an exponential function with the natural constant e as the base. The proportion of each exponential value to the sum of all exponential values ​​is then used as the contribution weight of each sequence boundary. The normalized exponential function can map any real value of significance score to the (0,1) interval while maintaining the relative order of the scores, ensuring that the sum of the contribution weights of all sequence boundaries is 1 and that each weight is positive, thus avoiding the subjectivity of human weighting. The significance score is the peak amplitude of the sequence boundary identification response value at the boundary location. The larger the peak amplitude, the more significant the geological characteristics of the sequence boundary, the stronger its control over the underlying strata, and the greater the contribution weight should be assigned.

[0024] Specifically, the expression for the three-dimensional spatial paleogeographic constraint field is: ; in, , , These represent the horizontal x-coordinate, horizontal y-coordinate, and continuous depth coordinates in three-dimensional space, respectively. This refers to the number of the layer sequence interface; This represents the total number of hierarchical interfaces. For three-dimensional paleogeographic constraint fields in coordinates The value at; For the first The depth constant value of each hierarchical interface; For the first The Gaussian decay width of each sequence interface, and ,in For the first The thickness difference between the strata above and below the sequence boundary. For all The largest stratigraphic thickness difference at each sequence boundary. The baseline Gaussian attenuation width, For adaptive adjustment index and For the first The contribution weight of each sequence interface is obtained by normalizing the index of the significance score of each sequence interface. For the natural constant An exponential function with base 0; coordinates The sedimentary facies type encoding value at the location; coordinates The thickness of the strata at that location; The auxiliary weighting coefficient for formation thickness; the continuous depth coordinate Z and the sequence boundary depth constant value. For depth parameters with the same dimensions, when the continuous depth coordinates The value is equal to the constant value of the interface depth of a certain sequence. At that time, the Gaussian attenuation term corresponding to the sequence boundary in the three-dimensional paleogeographic constraint field reaches its maximum value, indicating that the sequence boundary has the strongest control effect on the sedimentary environment at this depth; the Gaussian attenuation width is set. Difference in stratum thickness The positive correlation stems from the fact that the thickness difference between the strata above and below the sequence boundary reflects the magnitude of the sedimentary hiatus or sedimentary rate change corresponding to that boundary. A larger thickness difference indicates a more significant change in the sedimentary environment represented by that boundary, and its influence should be more widespread. Transforming stratigraphic thickness difference information into an adaptive adjustment of Gaussian attenuation width allows thicker interfaces to have a wider range of influence while thinner interfaces have a narrower range of influence. This reflects the difference in the sedimentary control effect of sequence interfaces at different levels. The value of is between 0 and 1, and is used to control the degree of nonlinear adjustment of the attenuation width by the thickness difference, so as to avoid the attenuation width being excessively amplified or reduced due to extreme thickness differences, and to ensure the robustness of the system.

[0025] The reference Gaussian attenuation width The method for determining this is as follows: The average thickness difference between the upper and lower strata above and below all sequence boundaries in the target area is statistically analyzed, and 0.5 times this average is taken as the threshold. The initial values ​​were then used, and sensitivity analysis was performed within the range of 0.1 to 1.0 times the average value through cross-validation. The value that resulted in the highest accuracy in identifying stratigraphic unit boundaries was selected as the final value. Value; the adaptive adjustment index Using grid search method The interval is scanned with a step size of 0.05, and the interval that minimizes the root mean square error of the prediction model on the validation set is selected. The above parameter determination method ensures that the values ​​of each parameter are based on data, avoiding the subjectivity of human experience in assigning values.

[0026] The weighting coefficients in the formula include two types: and , The contribution weight of each sequence interface, its subscript With sequence interface number Correspondingly, this indicates that different sequence interfaces have different contribution weights. In the summation term, each sequence interface is assigned a weight coefficient separately to reflect the different contributions of different sequence interfaces to paleogeographic constraints. These are auxiliary weighting coefficients for stratigraphic thickness. The subscript "3" indicates that, within the entire paleogeographic constraint field framework, the first set of weighting coefficients represents the contribution weight of sequence boundaries. (Corresponding to the sequence interface constraint in the first summation), the second set of weight coefficients is implicit in the Gaussian decay function of the sequence interface (corresponding to the width parameter of the Gaussian kernel). The third set of weighting coefficients are the auxiliary weighting coefficients for formation thickness. (Corresponding to the third auxiliary constraint), therefore the subscript "3" is used to indicate the position of the formation thickness in the overall constraint framework; this numbering system is only used to distinguish the weight coefficients of different sources; The value range of is [0,1], and the specific value is determined through cross-validation. When the stratigraphic thickness does not participate in the construction of the paleogeographic constraint field, The time indicates that the stratigraphic thickness and sequence boundary constraints have the same weight. In practical applications, it is usually taken to be between 0.2 and 0.5, so that the stratigraphic thickness participates in the calculation of the paleogeographic constraint field as an auxiliary correction term rather than a dominant term.

[0027] The stratigraphic analysis module further transforms stratigraphic boundary identification into solving for the local maximum surface of the three-dimensional paleogeographic constraint field along the depth direction. Specifically, it uses the depth points at various horizontal positions in three-dimensional space where the constraint field reaches a local maximum along the depth direction to construct the stratigraphic unit boundary surface, and delineates the top and bottom depths of each stratigraphic unit based on this surface. The principle behind transforming stratigraphic boundary identification into solving for the local maximum surface is that sequence boundaries are locations of abrupt changes in sedimentary environment, manifesting as local maximum points along the depth direction in the three-dimensional paleogeographic constraint field. This is because the sedimentary facies type encoding value changes abruptly at this location, and the Gaussian decay term of the sequence boundary reaches its peak. By solving for the local maximum surface of the constraint field, all stratigraphic boundary locations can be automatically identified without relying on manual interpretation of seismic profiles or well logging curves. This solution method does not require pre-setting the number of strata or stratification schemes; it is entirely data-driven and can automatically adapt to changes in geological characteristics of different regions, avoiding the inconsistencies in delineation caused by human factors in traditional methods.

[0028] Specifically, the stratigraphic analysis module includes a sequence boundary identification unit, a sedimentary facies analysis unit, and a paleogeographic pattern restoration unit. The sequence boundary identification unit calculates the difference between the natural gamma logging value and the mean of the sliding window, and the difference between the density logging value and the mean of the sliding window, at each depth point based on the well logging data and the seismic reflection data. The two values ​​are multiplied and then weighted within the sliding window to obtain a weighted average as the sequence boundary identification response value. The local peak position of the sequence boundary identification response value is taken as the sequence boundary position. The natural gamma logging value reflects the change in clay content in the formation, and the density logging value reflects the change in rock density. Sequence boundaries are usually accompanied by abrupt changes in clay content and rock density. The product of their differences can amplify the formation change signal and suppress noise. Weighted averaging smooths local fluctuations, making the response curve more stable. The sedimentary facies analysis unit uses... Based on the borehole lithology data and identified sequence boundaries, sedimentary facies analysis is performed on each sequence unit to generate sedimentary facies type coding values. The paleogeographic pattern restoration unit is used to perform Gaussian attenuation based on the sequence boundary identification results, sedimentary facies type coding values, and stratigraphic thickness, with the depth of each sequence boundary as the center and the Gaussian attenuation width determined by the difference in stratigraphic thickness above and below each sequence boundary. The attenuated sedimentary facies type coding values ​​are then weighted and superimposed with the contribution weight of each sequence boundary, and then weighted and summed with stratigraphic thickness as an auxiliary constraint term to establish the three-dimensional spatial paleogeographic constraint field. The basis for using stratigraphic thickness as an auxiliary constraint term in the weighted summation is that in areas with larger stratigraphic thickness, sedimentation is continuous and stable, and the paleogeographic environment has a more sufficient influence on stratigraphic development. Therefore, using stratigraphic thickness as a supplementary correction term of the paleogeographic constraint field can improve the accuracy of the constraint field in characterizing the sedimentary environment.

[0029] For sedimentary phosphate deposits in the target area, the mapping relationship of the sedimentary facies type coding values ​​is as follows: intertidal facies is coded as 5, subtidal facies as 4, platform facies as 3, slope facies as 2, and basin facies as 1. In this mapping relationship, the sedimentary facies type with a higher average grade has a larger coding value. Specifically, the coding method is as follows: Collect core descriptions and corresponding phosphate grade analysis data for all known phosphate deposits in the target area. After grouping by sedimentary facies type, calculate the arithmetic mean of the grades of all known phosphate deposits in each type. Sort the average grades from highest to lowest, assigning the sedimentary facies type with the highest average grade the largest coding value, the sedimentary facies type with the lowest average grade the smallest coding value, and assigning decreasing coding values ​​to intermediate types according to their average grades. For example, in a specific target area, if the average grade of the intertidal facies is statistically obtained as... The average grade of the intertidal facies is 18.5%, the average grade of the subtidal facies is 12.3%, the average grade of the platform facies is 8.7%, the average grade of the slope facies is 4.2%, and the average grade of the basin facies is 1.5%. Therefore, the average grades are assigned codes of 5, 4, 3, 2, and 1 respectively, from highest to lowest. According to the principles of phosphate sedimentology, the intertidal facies has high hydrodynamic conditions and frequent wave action, which is conducive to the enrichment and sorting of phosphorus particles, making it the most favorable sedimentary facies zone for phosphate deposits. The subtidal facies has relatively weak hydrodynamic conditions, but still possesses conditions for phosphorus enrichment. The platform facies has a stable sedimentary environment, and the phosphorus source is relatively limited. The slope and basin facies are far from the source area, resulting in insufficient phosphorus supply. Therefore, the intertidal facies has the highest code value, the basin facies has the lowest code value, and the intermediate types decrease in favor of mineralization. This coding order directly embeds the geological mineralization regularity into the numerical calculation framework.

[0030] The predictive modeling module is used to fit model parameters based on the known ore-controlling parameters and grade observations of phosphate deposits, extract the ore-controlling parameters of deep phosphate deposits in the target area, and generate three-dimensional distribution prediction results of phosphate resources through a nonlinear prediction model that includes power-law scaling, fracture index modulation, and direct coupling of three-dimensional paleogeographic constraint fields. The first nonlinear method is power-law scaling, which applies a power-law operation to the normalized value of each ore-controlling parameter to adjust its sensitivity. When the power-law is greater than 1, it amplifies the high-value area of ​​the ore-controlling parameter, while when it is less than 1, it compresses the low-value area. The basis for power-law scaling is that there is a significant threshold effect in the geological ore-controlling process. When the value of the ore-controlling parameter is below a certain critical value, its contribution to mineralization is weak and changes slowly. When the value of the ore-controlling parameter exceeds the critical value, its contribution to mineralization increases sharply. For example, after the thickness of the ore-bearing rock series reaches a certain value, the probability of mineralization does not increase nonlinearly but jumps.

[0031] Power exponent calculation can adaptively characterize this nonlinear response pattern of "gradual changes in low-value regions and dramatic changes in high-value regions." When the power exponent is less than 1, the curve is convex upwards, making it more sensitive to changes in low-value regions; when the power exponent is greater than 1, the curve is concave downwards, making it more sensitive to changes in high-value regions. The model training process automatically determines the optimal power exponent for each ore-controlling element, ensuring that it reflects the threshold effect while maintaining sufficient flexibility. The second nonlinear method is fracture exponent modulation, which uses the natural constant... The contribution of each ore-controlling parameter is modulated by the fracture density exponential function value. The basis for fracture exponential modulation is that the control effect of fracture structures on mineralization is not linearly cumulative, but exponentially cumulative. With each additional fracture, the connectivity and complexity of the fluid channel network expand exponentially, and the efficiency of ore-forming fluid transport increases exponentially. Using an exponential function to characterize this network effect is more in line with the principles of mineralization dynamics than linear superposition. The third nonlinear method is the direct coupling method of the three-dimensional paleogeographic constraint field, that is, the three-dimensional paleogeographic constraint field is directly superimposed into the prediction result, rather than being used as a normal input feature and combined with other ore-controlling parameters. The direct coupling method, which involves nonlinear transformation before prediction, is based on the following: Paleogeography is the primary geological factor controlling the formation of phosphate deposits. The distribution and grade of ore-bearing rock systems are fundamentally controlled by the paleogeographic pattern. If the paleogeographic constraint field and other ore-controlling elements undergo nonlinear transformations such as power-law scaling and fault-law modulation, their physical meaning and original spatial distribution characteristics may be distorted or weakened. Direct superposition can ensure that the paleogeographic constraint information participates in the prediction calculation in the form of the original numerical field throughout the process, avoiding the paleogeographic information being submerged by the weights of other ore-controlling elements in subsequent nonlinear transformations, thereby ensuring that the paleogeographic analysis results have a direct impact on the final prediction results.

[0032] The delineation module is used to establish an anisotropic three-dimensional variogram model and calculate the predicted value of the target location using Kriging interpolation based on this model. This predicted value is compared with the industrial grade index value for phosphate rock. Spatial locations with values ​​greater than or equal to the industrial grade index value are marked as predicted phosphate ore body boundaries, and a closed three-dimensional ore body envelope surface is generated. The industrial grade index value for phosphate rock refers to the minimum mineable grade in phosphate mining. This value is determined based on industrial indicators, and different types of phosphate rock have different industrial grade index values. For example, the boundary grade of phosphate rock type phosphate rock is usually [missing value]. content Industrial grade is usually The content is 15%~18%, and this index value can be preset in the system according to the specific mineral type and beneficiation conditions of the target area.

[0033] This embodiment also provides a three-dimensional prediction method for deep phosphate rock tectonic strata based on lithofacies paleogeographic constraints. The method includes the following steps: obtaining multi-source geological data for the target area, including borehole lithology data, well logging data, seismic reflection data, gravity data, and magnetic data; identifying sequence boundaries based on the multi-source geological data; performing Gaussian attenuation with the depth of each sequence boundary as the center and the Gaussian attenuation width determined by the difference in strata thickness above and below each sequence boundary; weighting and superimposing the attenuated sedimentary facies type codes with the contribution weights of each sequence boundary; and then weighting and summing the summations with strata thickness as an auxiliary constraint term to establish a three-dimensional spatial paleogeographic constraint field; and using the local maxima of this constraint field along the depth direction to construct the stratigraphic unit boundary surface. Based on this, the top and bottom depths of each stratigraphic unit are determined; the model parameters are fitted according to the known ore-controlling parameters and grade observation values ​​of phosphate deposits, and the normalized values ​​of each ore-controlling parameter are used as the base, and the power-law scaling is performed using the power-law adjustment coefficient as the exponent. The fracture index is modulated using the exponential function value of the fracture structure density, and the three-dimensional spatial paleogeographic constraint field is directly superimposed into the prediction results to establish a three-dimensional prediction model and generate a three-dimensional distribution prediction result of phosphate resources; an anisotropic three-dimensional variogram model is established, and the predicted value of the target location is calculated by Kriging interpolation. The predicted value is compared with the industrial grade index value of phosphate ore, and the spatial location that is greater than or equal to the index value is marked as the boundary of the predicted phosphate ore body, and a closed three-dimensional ore body envelope is generated.

[0034] Example 2: Based on Example 1, the expression for the three-dimensional distribution prediction result of the phosphate rock resources is as follows: ; in, , , These represent the horizontal x-coordinate, horizontal y-coordinate, and depth coordinates in three-dimensional space, respectively. coordinates Three-dimensional distribution prediction of phosphate resources at the location; Numbering of ore-controlling parameters; The total number of ore-controlling parameters; For the first Individual mineral control parameters in coordinates The normalized value at that point; coordinates Fault structure density at the location; For the first The power-law adjustment coefficient of each ore-controlling element parameter has an amplifying effect on the high-value region of the ore-controlling element when its value is greater than 1, and a compressing effect on the low-value region when its value is less than 1; exp(·) is an exponential function with the natural constant e as the base; For the first The fracture index adjustment coefficient of individual mineral control parameters; For the first The fusion weight of individual mineral control element parameters; The coupling coefficient of the three-dimensional paleogeographic constraint field is denoted as . The three-dimensional paleogeographic constraint field in coordinates ( The value at ().

[0035] The coupling coefficient satisfy ,in For the first The correlation coefficient between a ore-controlling element parameter and the three-dimensional paleogeographic constraint field; the meaning of the coupling coefficient calculation formula is: when the correlation coefficient between a ore-controlling element parameter and the paleogeographic constraint field is higher and its fusion weight is greater, the direct coupling strength of the paleogeographic constraint field also increases accordingly, so that the contribution of the paleogeographic constraint field and the overall contribution of each ore-controlling element remain consistent.

[0036] The power index adjustment coefficient The range of values ​​is ,in Take 0.1, Take 3.0, and The value range spans 1, allowing the model to learn both amplification and compression effects during training, ensuring a flexible nonlinear response capability to different ore-controlling factors; the fracture index adjustment coefficient The value of is non-negative to ensure that the predicted contribution increases monotonically with the increase of fault density, avoiding the physical paradox that the contribution decreases as the fault density increases, which is consistent with the geological law that faults promote mineralization; the fusion weight satisfy And each All are greater than or equal to 0. This constraint makes the contributions of each ore-controlling element comparable in a normalized way, preventing the weight of a certain element from being too large and dominating the prediction results. The above constraints encode geological knowledge into hard constraints of mathematical optimization problems, so that the optimization solution results are always within the geologically reasonable space.

[0037] The ore-controlling parameters extracted by the predictive modeling module include the thickness of the ore-bearing strata, the thickness of the phosphate rock layer, the grade of the phosphate rock layer, the dip angle of the strata, the density of fault structures, and the distance to paleogeographic facies zones. The thickness of the ore-bearing strata reflects the development degree of the phosphate-bearing strata; the greater the thickness, the longer the sedimentary process and the more sufficient the accumulation of ore-forming materials. The thickness of the phosphate rock layer directly reflects the scale of the ore body; the greater the thickness, the greater the resource quantity. The grade of the phosphate rock layer directly reflects the quality of the ore; the higher the grade, the greater the economic value. The dip angle of the strata reflects the degree of tectonic deformation; strata with a moderate dip angle are conducive to the migration and metasomatism of ore-forming fluids. The density of fault structures reflects the degree of fault development; the higher the density, the stronger the tectonic activity and the more developed the ore-forming fluid channels. The distance to paleogeographic facies zones reflects the distance of the target location from favorable sedimentary facies zones; the closer the distance, the greater the influence of favorable sedimentary environments. The above six parameters comprehensively characterize the phosphate mineralization conditions from six dimensions: sedimentary development, ore body scale, ore quality, tectonic deformation, fluid channels, and sedimentary environment.

[0038] The fracture structure density parameter is determined by three-dimensional Gaussian kernel fracture density calculation. Each fracture is considered a three-dimensional Gaussian distribution source. A three-dimensional Gaussian kernel is established centered on the fracture spatial coordinates. All Gaussian kernels are superimposed and divided by the kernel volume to obtain a three-dimensional continuous distribution of fracture structure density. Compared to the traditional method of summing fracture lengths per unit area, the three-dimensional Gaussian kernel fracture density calculation has the advantage of obtaining a smooth and continuous density field in three-dimensional space, avoiding abrupt density changes caused by fractures passing through grid boundaries. This makes the fracture structure density field differentiable in three-dimensional space, providing a mathematical basis for subsequent gradient-based optimization calculations. The paleogeographic facies distance parameter is determined by weighted facies distance calculation. For each sedimentary facies type, the shortest planar distance from the target location to the boundary of that type of facies is calculated, and then weighted by the ore-controlling weight coefficient of each facies type. This ore-controlling weight coefficient reflects the strength of different types of sedimentary facies in controlling mineralization; the closer the location is to a favorable facies, the greater the weight. Therefore, the influence of paleogeographic facies is reflected in the prediction results through distance attenuation.

[0039] The predictive modeling module also includes a model iteration unit, which iteratively optimizes the joint optimization objective function using the alternating direction multiplier method. Each iteration sequentially updates the fusion weights, power exponent adjustment coefficients, and fracture exponent adjustment coefficients. The fusion weight update employs the proximal gradient operator to ensure its non-negativity and the constraint that its summation equals 1. The reason for using the alternating direction multiplier method is that the joint optimization objective function contains three regularization terms and multiple constraints, belonging to a constrained non-convex optimization problem. The alternating direction multiplier method can decompose the original problem into three sub-problems that are solved separately, with each iteration updating only one type of parameter while fixing the other two types. This reduces the difficulty of solving the problem and ensures convergence. In each iteration, the power exponent adjustment coefficient and the break exponent adjustment coefficient are first fixed, and the fusion weights are updated to ensure that they satisfy the constraints of non-negativity and summation to 1. The proximal gradient operator achieves constraint satisfaction by projecting onto the feasible region. Then, the fusion weights and the break exponent adjustment coefficient are fixed, and the power exponent adjustment coefficient is updated to ensure that it falls within a preset interval. Finally, the fusion weights and the power exponent adjustment coefficient are fixed, and the break exponent adjustment coefficient is updated to ensure that it is non-negative. This process is repeated iteratively until convergence. The convergence criterion is that the infinite norm of the three types of parameter vectors in two consecutive iterations is less than a preset threshold. .

[0040] The expression for the anisotropic three-dimensional variogram model is as follows: ; in, It is a three-dimensional spatial lag distance vector. These represent the components of the vector in the horizontal x-axis direction, the horizontal y-axis direction, and the depth direction, respectively. The coordinates mentioned above are in the dimension of depth. The hysteresis distance vector The corresponding variogram value; This is the base value; For variable range; It is an anisotropic tensor matrix; For the natural constant An exponential function with base 0; The hysteresis distance vector Transpose of; It is a quadratic form of the hysteresis distance vector and the anisotropic tensor matrix; the anisotropic tensor matrix The orthogonal rotation matrix is ​​obtained by eigenvalue decomposition of the diagonal eigenvalue matrix and the orthogonal rotation matrix, wherein the orthogonal rotation matrix is ​​the gradient vector field formed by taking the partial derivatives of the three-dimensional paleogeographic constraint field along each coordinate axis. Principal component analysis (PCA) is used to determine the gradient vector, specifically by using its three components as three feature dimensions. A data matrix is ​​constructed from the three component values ​​of all grid points, and PCA is performed to obtain the three principal component directions as the three column vectors of the orthogonal rotation matrix. The three diagonal elements of the diagonal eigenvalue matrix are the reciprocals of the squares of the range parameters of the three principal directions. The gradient vector field... The expression is The basis for using anisotropy lies in the fact that phosphate deposits are controlled by paleocurrent direction and basin morphology, exhibiting significantly different spatial correlations along and perpendicular to the paleocurrent direction. The anisotropic tensor matrix can accurately characterize the directional spatial structure of phosphate grade. The principle of determining the orthogonal rotation matrix through principal component analysis is that the gradient direction field of the three-dimensional paleogeographic constraint field reflects the spatial change direction of the paleogeographic pattern. The direction with the largest gradient change corresponds to the direction of the most drastic change in the paleogeographic pattern, while the direction with the smallest gradient change corresponds to the direction of the most gradual change in the paleogeographic pattern. The direction with the largest variance extracted by principal component analysis is the principal direction of the anisotropic tensor, and the range parameter reflects the spatial correlation range along this direction. The above-mentioned direct conversion of the directional information of the paleogeographic pattern into the anisotropic parameters of the variogram model realizes automated anisotropic interpolation driven by geological constraints, eliminating the need for geologists to manually specify the anisotropic direction based on experience, and avoiding the influence of subjective factors on the boundary delineation results. (Platform value) The range was determined by fitting the variogram function experimentally. It is obtained by least-squares fitting of the range parameters in the variogram curve.

[0041] When the delineation module calculates the predicted value of the target location using Kriging interpolation, the Kriging interpolation weight coefficients are determined by solving the Kriging equations. These equations include the variogram values ​​between each known phosphate rock point and all other known phosphate rock points, the variogram values ​​between each known phosphate rock point and the target location, and the constraint that the weight coefficients sum to 1. By solving these equations, the optimal weight coefficients for each known phosphate rock point are obtained, thus yielding an unbiased optimal estimate of the target location. The advantages of Kriging interpolation compared to traditional interpolation methods such as inverse distance weighting are: Kriging interpolation characterizes the spatial structure of the data through the variogram, and the determination of the weight coefficients is based on the principle of statistical optimality, ensuring that the interpolation result is unbiased in mathematical expectation and has the smallest variance, while also providing an estimate of the uncertainty of the interpolation result; coupling Kriging interpolation with anisotropic variograms allows the boundary delineation to simultaneously possess geological directional constraints and statistical optimality.

[0042] The system also includes a 3D display module and a prediction verification module. The 3D display module is used to generate a 3D geological model and a 3D phosphate resource distribution prediction map of the target area based on the 3D ore body envelope. The 3D geological model displays the spatial distribution and structural morphology of stratigraphic units in a 3D form, and the 3D phosphate resource distribution prediction map displays the distribution of phosphate resource prediction index values ​​in 3D space in the form of color gradients or isosurfaces, providing users with a visualization tool to intuitively judge favorable mineralization areas. The prediction verification module is used to determine the preferred exploration target area based on the 3D distribution prediction results of phosphate resources and the 3D ore body envelope, and outputs the 3D spatial coordinate range of the preferred exploration target area. The method for determining the preferred exploration target area is to select several closed areas with the highest prediction values ​​as the preferred target area. The center coordinates and boundary range of the closed area are the 3D spatial coordinate range. This output can be directly used to guide the deployment of subsequent drilling projects.

[0043] This embodiment also provides a three-dimensional prediction method for deep phosphate rock tectonic stratigraphy based on lithofacies paleogeographic constraints. The method includes the following steps: obtaining multi-source geological data for the target area, including borehole lithology data, well logging data, seismic reflection data, gravity data, and magnetic data; identifying sequence boundaries based on the multi-source geological data; performing Gaussian attenuation with the depth of each sequence boundary as the center and the Gaussian attenuation width determined by the difference in strata thickness above and below each sequence boundary; weighting and superimposing the attenuated sedimentary facies type codes with the contribution weights of each sequence boundary; and then weighting and summing the results with strata thickness as an auxiliary constraint term to establish a three-dimensional spatial paleogeographic constraint field; using the local maxima of this constraint field along the depth direction to construct the stratigraphic unit boundary surface, and delineating the top and bottom depths of each stratigraphic unit accordingly; and determining the ore-controlling factors of known phosphate deposits. The model parameters are fitted with the element parameters and grade observation values. The normalized values ​​of each ore-controlling element parameter are used as the base, and a power-law scaling factor is used as the exponent. The fault index is modulated using the exponential function value of the fault structure density. The three-dimensional paleogeographic constraint field is directly superimposed onto the prediction results to establish a three-dimensional prediction model and generate a three-dimensional distribution prediction result for phosphate resources. An anisotropic three-dimensional variogram model is established, and the predicted value of the target location is calculated using Kriging interpolation. This predicted value is compared with the phosphate industrial grade index value. Spatial locations greater than or equal to this index value are marked as the predicted phosphate ore body boundary, and a closed three-dimensional ore body envelope is generated. The specific construction methods of the three-dimensional paleogeographic constraint field, the three-dimensional prediction model, and the anisotropic three-dimensional variogram model in the above steps are the same as those described in Examples 1 and 2.

[0044] Example 3: Based on Examples 1 and 2, a sedimentary phosphate mining exploration area was selected as the target area. The known borehole data covers an area of ​​approximately 35 square kilometers, with a total of 76 boreholes. Among them, 32 phosphate mining points are known. Geological surveys have revealed that the area has a well-developed phosphate-bearing rock series. The phosphate layer is hosted in the middle section of the Doushantuo Formation of the Sinian System. The main sedimentary facies types include five types: intertidal facies, subtidal facies, platform facies, slope facies, and basin facies.

[0045] Following the method described in Example 1, multi-source geological data for the area were first obtained, including lithological logging data and well logging curve data from 76 boreholes, three-dimensional seismic reflection data covering the entire area, and gravity and magnetic measurement data at a scale of 1:50,000. Sequence boundaries were identified based on the well logging curve data, and a total of four main sequence boundaries were identified. =4), depth values ​​of each sequence interface to The corresponding thickness differences between the upper and lower strata are 320m, 480m, 650m, and 820m, respectively. to The m lengths are 45m, 62m, 38m, and 53m respectively. It is 62m; the reference Gaussian attenuation width is taken. Adaptive adjustment index The Gaussian decay width of each sequence interface was calculated. to The depths are 26.3m, 30.0m, 24.1m, and 28.2m, respectively. Based on the mapping relationship of sedimentary facies type coding values, the average grades of 32 known phosphate deposits were statistically analyzed: the average grade for intertidal facies was 18.5%, for subtidal facies it was 12.3%, for platform facies it was 8.7%, for slope facies it was 4.2%, and for basin facies it was 1.5%, assigned coding values ​​of 5, 4, 3, 2, and 1, respectively. The contribution weights of each sequence boundary were then used to determine the appropriate grades. The attenuated sedimentary facies type codes are weighted and superimposed, and then supplemented with the stratigraphic thickness auxiliary term ( We take a weighted sum of 0.3 and establish a three-dimensional paleogeographic constraint field.

[0046] The local maximum surface of the three-dimensional paleogeographic constraint field along the depth direction is solved, and four stratigraphic unit boundary surfaces are automatically identified. Based on these, the top and bottom depths of five stratigraphic units are delineated. Compared with the results of manual interpretation, the average absolute error of stratigraphic boundary identification is 6.8m, and the maximum error is 12.5m, indicating that the automatic identification method of the present invention has high accuracy.

[0047] Following the method described in Example 2, six types of ore-controlling parameters were extracted. The 32 known phosphate deposits were divided into a training set (22 deposits) and a validation set (10 deposits) at a 7:3 ratio. The joint objective function was iteratively optimized using the alternating direction multiplier method, converging after 156 iterations. The power-law adjustment coefficient was then used. The optimization results were 1.32, 0.78, 1.54, 0.95, 2.01, and 1.17, respectively, and the fracture index adjustment coefficient was... The optimization results were 0.21, 0.35, 0.18, 0.42, 0.56, and 0.29, respectively, with fusion weights... The optimization results are as follows: and 0.10; coupling coefficient according to The value is calculated to be 0.35; a nonlinear prediction model is established to generate the three-dimensional distribution prediction results of phosphate rock resources.

[0048] An anisotropic three-dimensional variogram model was established. Principal component analysis determined the range parameters for the three principal directions to be 240m, 180m, and 120m, respectively, and the sill value was... The fit is 0.82, and the range is [not specified]. The fitted value is 150m; the predicted value of the target location is calculated by Kriging interpolation, and the industrial grade index of phosphate rock is taken as 150m. With a content of 12%, spatial locations where the predicted value is greater than or equal to the predicted value are marked as the boundaries of the predicted phosphate ore body, generating a closed three-dimensional ore body envelope.

[0049] In practical application, based on existing technology, the three-dimensional prediction of deep phosphate rock strata is divided into four sequentially connected stages: data preparation stage, paleogeographic constraint field construction and stratigraphic division stage, nonlinear prediction modeling stage, and three-dimensional ore body delineation and target area output stage. Each stage is executed sequentially, and the output of the previous stage serves as the input of the next stage.

[0050] Data preparation phase: Personnel input the acquired borehole lithology data, well logging data, seismic reflection data, gravity data, and magnetic data for the target area into the data input module. For borehole lithology data, lithology names and codes are standardized according to stratigraphic logging specifications. For well logging data, depth and environmental corrections are performed on natural gamma curves and density curves to eliminate the influence of wellbore variation and mud intrusion on the measured values. For seismic reflection data, time-depth conversion is performed to convert the time-domain reflection profile into a depth-domain reflection profile. For gravity and magnetic data, topographic and diurnal variation corrections are performed to obtain Bouguer gravity anomaly data and magnetic susceptibility data. After preprocessing, all data is imported into the data input module. If core recovery rate data for some boreholes is missing in the target area, the average core recovery rate of adjacent boreholes can be used temporarily and updated as new data becomes available.

[0051] In the paleogeographic constraint field construction and stratigraphic division stage: The implementers specify the sliding window width parameter required for sequence boundary identification. This parameter is set based on the average stratigraphic thickness of the target area, generally ranging from 10m to 30m. If the target area's strata are predominantly thin (single layer thickness less than 5m), a smaller window width (10m to 15m) is used; if predominantly thick (single layer thickness greater than 10m), a larger window width (20m to 30m) is used. After the stratigraphic analysis module is activated, the sequence boundary identification unit automatically processes the well logging data, calculating the difference between the natural gamma logging value and the window mean, and the difference between the density logging value and the window mean at each depth point. After multiplication, a weighted average is taken within the window to obtain the sequence boundary identification response curve for the entire well section. The local peak positions of the response curve are automatically marked as sequence boundary locations. Subsequently, the sedimentary facies analysis unit determines the sedimentary facies of each stratigraphic unit based on borehole lithology data. The implementation personnel verify and correct the automatically determined results to ensure the accuracy of sedimentary facies type attribution. The paleogeographic pattern restoration unit automatically establishes a three-dimensional spatial paleogeographic constraint field based on the sequence boundary identification results and sedimentary facies type coding values. During the calculation process, the implementation personnel can adjust the values ​​according to the regional geological understanding. If the stratigraphic thickness in the target area varies significantly (thickness variation coefficient greater than 0.3), the values ​​can be appropriately increased. The thickness should be increased to 0.4-0.5 to enhance the correction effect of stratigraphic thickness on the paleogeographic constraint field. If the stratigraphic thickness changes gradually (thickness variation coefficient less than 0.15), it can be appropriately reduced. The depth is set to 0.1~0.2. After the three-dimensional paleogeographic constraint field is established, the stratigraphic analysis module automatically searches for local maxima along the depth direction to form the boundary surface of the stratigraphic unit, and delineates the top and bottom depths of each stratigraphic unit. The implementers review the automatically delineated stratigraphic boundaries and manually adjust any unreasonable boundaries. After adjustment, the module automatically updates the thickness parameters of adjacent stratigraphic units.

[0052] In the nonlinear prediction modeling phase: The implementers specify the ratio of training and validation sets for known phosphate deposits. A 7:3 ratio is used when there are more than 30 known phosphate deposits; an 8:2 ratio is used when there are 15 to 30; and leave-one-out cross-validation is used when there are fewer than 15. The prediction modeling module automatically extracts the spatial distribution data of each parameter based on the ore-controlling element parameter types specified by the implementers (selecting all or some types from six categories: ore-bearing rock thickness, phosphate layer thickness, phosphate layer grade, stratigraphic dip angle, fault density, and paleogeographic facies distance). Subsequently, the model iteration unit initiates an automatic optimization process, using the alternating direction multiplier method to jointly optimize the three types of parameters: fusion weights, power exponent adjustment coefficients, and fault exponent adjustment coefficients. During the optimization process, [further details are needed]. Personnel monitor the loss function's decline curve in real time to ensure normal convergence of the iterations. If the loss function decreases by less than 0.1% of the initial loss value in 20 consecutive iterations, the iterations are considered converged. After optimization, personnel verify the prediction accuracy on the validation set. If the validation set accuracy meets the requirements (root mean square error less than 0.15 or coefficient of determination greater than 0.75), the model parameters are confirmed and locked. If the accuracy does not meet the requirements, the regularization coefficient is adjusted and retraining is performed, or the ratio of the training set to the validation set is adjusted and retraining is performed. After model confirmation, the prediction modeling module generates a three-dimensional spatial prediction result for the entire target area. This result is stored in the form of a three-dimensional grid, with each grid node corresponding to a prediction index value. The relative magnitude of the prediction index value reflects the relative probability of phosphate mine development.

[0053] 3D Ore Body Delineation and Target Area Output Stage: Based on the industrial grade standards of the target mining area, the implementers input the industrial grade index value of phosphate rock through the delineation module. This value is determined according to the "Specifications for Geological Exploration of Phosphate Rocks" (DZ / T 0209-2020) and the deposit type. For example, the boundary grade of a phosphate rock-type phosphate rock deposit is... Content 8%~10%, industrial grade The content is 15%~18%; the delineation module automatically calculates the Kriging interpolation prediction value of each grid node in three-dimensional space, marks the continuous area with the prediction value greater than or equal to the industrial grade index value as the predicted ore body, and automatically tracks the boundary nodes of the predicted ore body. The moving cube algorithm is used to generate a closed three-dimensional ore body envelope surface; the three-dimensional display module overlays the ore body envelope surface with the three-dimensional geological model, and implementers can rotate, zoom, and cut in the three-dimensional view to observe the spatial morphology of the ore body and its spatial relationship with the strata; the prediction verification module automatically identifies several closed areas with the highest prediction index values ​​as preferred exploration target areas, and outputs the center coordinates and boundary range of each target area in descending order of prediction index values. Implementers formulate drilling verification plans based on the target area ranking results, and prioritize the verification project for the target area with the highest prediction index value.

[0054] The above description is a further detailed explanation of the invention in conjunction with specific embodiments, and it should not be considered that the specific embodiments of the invention are limited to this. For those skilled in the art to which this invention pertains and related fields, any extensions, operation methods, and data substitutions made based on the technical solution concept of this invention should fall within the protection scope of this invention.

Claims

1. A three-dimensional prediction system for the tectonic stratigraphy of deep phosphate deposits based on lithofacies paleogeographic constraints, characterized in that, include: The data input module obtains multi-source geological data for the target area; The stratigraphic analysis module identifies the sequence boundaries of the target area based on the multi-source geological data. It performs Gaussian attenuation with the depth of each sequence boundary as the center and the Gaussian attenuation width determined by the difference in strata thickness above and below each sequence boundary. The attenuated sedimentary facies type codes are then weighted and superimposed with the contribution weight of each sequence boundary, and finally weighted and summed with strata thickness as an auxiliary constraint term to establish a three-dimensional paleogeographic constraint field. The stratigraphic analysis module also transforms stratigraphic boundary identification into solving the local maximum surface of the three-dimensional paleogeographic constraint field along the depth direction. That is, the depth points at each horizontal position in the three-dimensional space where the three-dimensional paleogeographic constraint field reaches a local maximum along the depth direction are used to form the stratigraphic unit boundary surface, and the top and bottom depths of each stratigraphic unit are determined according to the stratigraphic unit boundary surface. The prediction modeling module fits model parameters based on known phosphate deposit control parameters and grade observation values, extracts phosphate control parameters for deep phosphate deposits in the target area, and generates three-dimensional distribution prediction results for phosphate resources through a nonlinear prediction model that includes power scaling, fracture index modulation, and direct coupling of three-dimensional paleogeographic constraint fields. The delineation module establishes an anisotropic three-dimensional variogram model. Using the predicted three-dimensional distribution of phosphate resources as known sampling point data, the module calculates the predicted Kriging interpolation value through Kriging interpolation. The predicted Kriging interpolation value is compared with the industrial grade index value of phosphate rock. Spatial locations that are greater than or equal to the industrial grade index value of phosphate rock are marked as the predicted phosphate ore body boundary, and a closed three-dimensional ore body envelope surface is generated.

2. The system according to claim 1, characterized in that, The expression for the three-dimensional paleogeographic constraint field is: ； in, , , These represent the horizontal x-coordinate, horizontal y-coordinate, and continuous depth coordinates in three-dimensional space, respectively. This refers to the number of the layer sequence interface; This represents the total number of hierarchical interfaces. For three-dimensional paleogeographic constraint fields in coordinates The value at; For the first The depth constant value of each layer interface; For the first The Gaussian decay width of each sequence interface, and . ,in For the first The thickness difference between the strata above and below the sequence boundary. For all The largest stratigraphic thickness difference at each sequence boundary. The baseline Gaussian attenuation width, For adaptive adjustment index and ; For the first The contribution weights of each sequence interface; exp(·) is a function of natural constants. An exponential function with base 0; coordinates ( The sedimentary facies type code value at () is used. coordinates The thickness of the strata at that location; This is an auxiliary weighting coefficient for formation thickness.

3. The system according to claim 1, characterized in that, The multi-source geological data includes borehole lithology data, well logging curve data, seismic reflection data, gravity data, and magnetic data; the stratigraphic analysis module includes: The sequence interface identification unit calculates the difference between the natural gamma logging value and the mean of the sliding window, and the difference between the density logging value and the mean of the sliding window at each depth point based on the well logging data and the seismic reflection data. The two values ​​are multiplied and the weighted average is taken within the sliding window as the sequence interface identification response value. The local peak position of the sequence interface identification response value is taken as the sequence interface position. The sedimentary facies analysis unit performs sedimentary facies analysis on each sequence unit based on the borehole lithology data and the identified sequence boundaries, and generates sedimentary facies type coding values. The paleogeographic pattern restoration unit, based on the sequence boundary identification results, sedimentary facies type coding values, and stratigraphic thickness, performs Gaussian attenuation with the depth of each sequence boundary as the center and the Gaussian attenuation width determined by the difference in stratigraphic thickness above and below each sequence boundary. The attenuated sedimentary facies type coding values ​​are then weighted and superimposed with the contribution weight of each sequence boundary, and finally weighted and summed with stratigraphic thickness as an auxiliary constraint term to establish the three-dimensional spatial paleogeographic constraint field.

4. The system according to claim 1, characterized in that, For sedimentary phosphate deposits in the target area, the mapping relationship of the sedimentary facies type coding values ​​is as follows: intertidal facies is coded as 5, subtidal facies as 4, platform facies as 3, slope facies as 2, and basin facies as 1.

5. The system according to claim 2, characterized in that, The expression for the three-dimensional distribution prediction result of the phosphate rock resources is as follows: ； in, , , These represent the horizontal x-coordinate, horizontal y-coordinate, and depth coordinates in three-dimensional space, respectively. coordinates Three-dimensional distribution prediction of phosphate rock resources at the location; Numbering of ore-controlling parameters; The total number of ore-controlling parameters; For the first Individual mineral control parameters in coordinates The normalized value at that point; coordinates Fault structure density at the location; For the first The power-law adjustment coefficients of the mineral control parameters; exp(·) is an exponential function with the natural constant e as the base; For the first The fracture index adjustment coefficient of individual mineral control parameters; For the first The fusion weight of individual mineral control element parameters; The coupling coefficient of the three-dimensional spatial paleogeographic constraint field; The three-dimensional paleogeographic constraint field in coordinates The value at the location; the coupling coefficient satisfy ,in For the first Correlation coefficients between individual mineral control parameters and the three-dimensional paleogeographic constraint field; The power index adjustment coefficient The range of values ​​is and The fracture index adjustment coefficient The value of is a non-negative number; the fusion weight satisfy And each All are greater than or equal to 0.

6. The system according to claim 5, characterized in that, The ore-controlling parameters extracted by the prediction modeling module include the thickness of the ore-bearing rock series, the thickness of the phosphate rock layer, the grade of the phosphate rock layer, the dip angle of the strata, the density of fault structures, and the distance of paleogeographic facies zones. The fracture structure density parameter is determined by the three-dimensional Gaussian kernel fracture density calculation, that is, each fracture is taken as a three-dimensional Gaussian distribution source, a three-dimensional Gaussian kernel is established with the fracture spatial coordinates as the center, and all Gaussian kernels are superimposed and divided by the Gaussian kernel volume to obtain the three-dimensional continuous distribution of fracture structure density. The paleogeographic facies distance parameter is determined by weighted facies distance calculation, that is, for each sedimentary facies type, the shortest planar distance from the target location to the boundary of that sedimentary facies type is calculated, and the distance is summed by weighting the ore-controlling weight coefficients of each facies type.

7. The system according to claim 5, characterized in that, The prediction modeling module also includes a model iteration unit, which iteratively optimizes the joint optimization objective function using the alternating direction multiplier method. In each iteration, the fusion weight, the power exponent adjustment coefficient, and the fracture exponent adjustment coefficient are updated sequentially. The fusion weight update uses the proximal gradient operator to ensure its non-negativity and the constraint that the summation is 1.

8. The system according to claim 2, characterized in that, The expression for the anisotropic three-dimensional variogram model is: ； in, It is a three-dimensional spatial lag distance vector. , , These represent the components of the lag distance vector in three-dimensional space along the horizontal x-axis, horizontal y-axis, and depth directions, respectively. It belongs to the same depth dimension as the depth coordinate Z; The hysteresis distance vector The corresponding variogram value; This is the base value; For variable range; It is an anisotropic tensor matrix; For the natural constant An exponential function with base 0; The hysteresis distance vector transpose; It is a quadratic form of the hysteresis distance vector and the anisotropic tensor matrix.

9. The system according to claim 1, characterized in that, Also includes: The 3D display module generates a 3D geological model and a 3D phosphate resource distribution prediction map of the target area based on the 3D ore body envelope. The prediction and verification module determines the preferred exploration target area based on the predicted three-dimensional distribution of the phosphate rock resources and the three-dimensional ore body envelope, and outputs the three-dimensional spatial coordinate range of the preferred exploration target area.

10. The prediction method for the system according to any one of claims 1 to 9, characterized in that, Includes the following steps: Step 1: Obtain multi-source geological data for the target area, including borehole lithology data, well logging curve data, seismic reflection data, gravity data, and magnetic data; Step 2: Identify sequence boundaries based on the multi-source geological data. Using the depth of each sequence boundary as the center and the Gaussian attenuation width determined by the difference in strata thickness above and below each sequence boundary, perform Gaussian attenuation. Then, weight the attenuated sedimentary facies type codes with the contribution weight of each sequence boundary and sum them up with strata thickness as an auxiliary constraint term to establish a three-dimensional paleogeographic constraint field. The local maxima of the three-dimensional paleogeographic constraint field along the depth direction constitute the stratigraphic unit boundary surface, and the top and bottom depths of each stratigraphic unit are delineated accordingly. Step 3: Fit model parameters based on known ore-controlling parameters and grade observation values ​​of phosphate deposits. Use the normalized values ​​of each ore-controlling parameter as the base and the power-law adjustment coefficient as the exponent for power-law scaling. Use the exponential function value of the fracture structure density for fracture-law modulation. Then, directly superimpose the three-dimensional paleogeographic constraint field into the prediction results to establish a three-dimensional prediction model and generate a three-dimensional distribution prediction result of phosphate resources. Step 4: Establish an anisotropic three-dimensional variogram model. Using the predicted three-dimensional distribution of phosphate resources as known sampling point data, calculate the predicted Kriging interpolation value through Kriging interpolation. Compare the predicted Kriging interpolation value with the industrial grade index value of phosphate rock. Mark the spatial locations that are greater than or equal to the industrial grade index value of phosphate rock as the predicted phosphate ore body boundary and generate a closed three-dimensional ore body envelope surface.

Images