Self-adaptive multi-modal feature fusion well logging interpretation method based on multi-task joint learning

Abstract

The invention provides an adaptive multi-modal feature fusion well logging interpretation method based on multi-task joint learning, and the model comprehensively employs a graph convolution network, a multi-head attention mechanism, a gating mechanism and an adaptive loss weight adjustment algorithm, and effectively integrates multi-well group well logging data and geological structure information. And the complex spatial relationship and geological similarity between the wells are fully captured. The model can process multi-source data in a unified manner, extracts high-level feature representation containing spatio-temporal correlation, avoids repeated feature learning, reduces calculation complexity, and ensures the consistency and complementarity of prediction task features at the same time. The self-adaptive loss weight adjustment mechanism dynamically optimizes weight distribution according to the loss change of each task, realizes balance optimization among different tasks, and comprehensively improves the prediction precision of each task. The physical constraint regularization module ensures that the prediction result not only has statistical accuracy, but also conforms to the actual geological physical law, so that the reasonability and reliability of the prediction result are improved.

Description

Technical Field

[0001] The present invention relates to the fields of geological engineering and oil and gas exploration, and more specifically to a well logging interpretation method based on multi-task joint learning and adaptive multi-modal feature fusion. Background Art

[0002] In the field of geological engineering and oil and gas exploration, accurate prediction of well permeability (PEM), porosity (POR), formation saturation (SOI) and lithology classification is of great significance for reservoir evaluation, geological modeling and resource development. Traditional logging data analysis methods mainly rely on single-task learning models, which usually predict a single indicator, such as permeability or porosity, but have difficulty in processing multiple related tasks at the same time. This single-task analysis method is not only inefficient, but also fails to fully utilize the inherent correlation between different logging indicators, resulting in limited comprehensiveness and accuracy of the prediction results. In addition, when processing complex geological data, single-task models often ignore the spatial relationship and geological similarity between wells, making it difficult to capture the impact of geological structure on logging indicators, thus affecting the overall prediction performance.

[0003] With the development of machine learning technology, multi-task learning (MTL) has gradually become an effective means to solve the above problems. Multi-task learning can share information between multiple related tasks by sharing feature extraction layers, thereby improving the generalization ability and prediction accuracy of the overall model. However, existing multi-task learning models still face many challenges when processing geological logging data. First, geological logging data has a high degree of spatial and depth dependence, and traditional multi-task learning models are difficult to fully capture the complex spatial relationships and geological similarities between wells. Secondly, there are significant differences in the loss scale and learning rate of different logging indicators. Directly adding the loss functions of each task may cause the loss of a certain task to dominate the overall optimization process, thereby affecting the prediction accuracy of other tasks. In addition, geological logging data itself is subject to physical laws and geological phenomena, such as porosity decreases with increasing depth, water saturation decreases with increasing depth, etc. Existing multi-task learning models often lack a mechanism to incorporate these physical constraints into the model training process, resulting in insufficient rationality and reliability of prediction results in practical applications.

[0004] Graph Convolutional Networks (GCN), as a deep learning model that can effectively process graph structured data, has achieved remarkable results in many fields in recent years. GCN can capture the spatial relationship and structural information between nodes by performing convolution operations on the graph structure, and is suitable for processing geographical and geological similarities between wells. However, relying solely on GCN for multi-task learning still faces the challenges of the above-mentioned multi-task learning models, especially in how to effectively integrate multiple adjacency matrix information and handle the loss balance between different tasks. At the same time, existing technologies such as multi-head attention mechanism and gating mechanism can better improve the model feature selectivity and expression ability, but in existing research, there is no case of fully utilizing these mechanisms to optimize multi-task learning models, especially the research on combined prediction with geophysical constraints is relatively small.

[0005] Geological logging data often face the problem of inconsistent and missing data formats, and data preprocessing is relatively complex. Traditional data filling methods, such as mean filling or interpolation, can solve the problem of missing data to a certain extent, but they often cannot fully retain the geological characteristics and spatial correlation of the data, affecting the training effect and prediction performance of subsequent models. The use of data enhancement technology can effectively improve the robustness of the model and adapt to sparse data, but the existing enhancement methods have not been fully developed in the application of geological logging data. In particular, how to perform multi-dimensional data enhancement while maintaining geological rationality is still a difficult problem that needs to be solved urgently. Summary of the invention

[0006] In response to the above problems, a new model that can handle multiple tasks simultaneously, effectively integrate multimodal features and take into account geophysical laws is urgently needed. This model not only needs to optimize the prediction accuracy of each task under the multi-task learning framework, but also needs to introduce a physical constraint regularization module to ensure that the prediction results are consistent with the actual geophysical laws on the basis of statistical accuracy. In addition, the model should have advanced data preprocessing and enhancement capabilities to solve the problems of inconsistent and missing data formats, and generate virtual well data that conforms to geological laws through multi-dimensional data enhancement methods to improve the model's robustness and adaptability to sparse data.

[0007] The present invention proposes an adaptive multimodal feature fusion model based on multi-task joint learning, which comprehensively uses graph convolutional networks, multi-head attention mechanisms, gating mechanisms and adaptive loss weight adjustment algorithms to achieve high-precision joint prediction of permeability, porosity, formation saturation and lithology classification. The model fully integrates logging data and geological structure information through a shared feature extraction layer, extracts high-level feature representations containing time and space correlations, avoids repeated learning of features for each task, and reduces computational complexity. At the same time, the model introduces a multi-task loss balance module, which adaptively adjusts the loss weights of each task to ensure that different tasks achieve balanced optimization during the training process, prevents a certain task from dominating the overall optimization process, and thus improves the prediction accuracy of all tasks. In addition, the physical constraint regularization module incorporates geophysical laws into the model training process, and by constraining geological relationships such as porosity and depth, water saturation and depth, porosity and density, and water saturation and relative permeability, ensures that the prediction results are in line with the actual geophysical laws on the basis of statistical accuracy, and improves the rationality and reliability of the prediction results.

[0008] The technical solution of the present invention is specifically as follows:

[0009] A well logging interpretation method based on adaptive multi-modal feature fusion of multi-task joint learning includes the following steps:

[0010] S1: Data input and preprocessing stage, including obtaining multi-well group logging data, well location data and geological structure data, and performing format conversion, missing value filling, coordinate conversion and spatial connection on the data to ensure data integrity and consistency;

[0011] S2: Data enhancement processing stage, using a multi-dimensional hybrid data enhancement method, including spatial enhancement, feature enhancement, and depth enhancement, to generate virtual well data that conforms to geological laws, improving the model's robustness and adaptability to sparse data;

[0012] S3: Dynamic adjacency matrix construction: Based on the geographic spatial location, geological characteristics similarity and depth sequence similarity of the wells, the spatial adjacency matrix, geological similarity adjacency matrix and depth sequence adjacency matrix are constructed respectively to accurately capture the complex spatial and geological relationships between wells;

[0013] S4: In the adaptive multimodal feature fusion model construction phase, graph convolution, multi-head attention mechanism and gating mechanism are used to extract high-level features from the fused features, and the information of various adjacency matrices is integrated through dynamic weighting.

[0014] Generate unified feature representation;

[0015] S106 multi-task joint learning achieves simultaneous high-precision prediction of PEM, POR, SOI and lithology classification through shared feature extraction layer, multi-task loss balance module and multi-task prediction module.

[0016] S5: Model training, using the preprocessed and enhanced feature matrix and the constructed adjacency matrix to train the constructed multimodal feature fusion model, adopting a multi-task joint learning framework, combining cross entropy loss and mean square error loss with physical constraint regularization terms to ensure that the prediction results conform to the laws of geophysics.

[0017] The present invention significantly improves the prediction accuracy and reliability of key logging indicators in oil and gas exploration through systematic implementation steps and combined with multiple innovative technologies. The main innovations and advantages of the present invention are summarized below in the order of implementation steps:

[0018] 1. Efficient data input and preprocessing

[0019] The present invention proposes a complete data preprocessing method to solve the problem of inconsistent and missing logging data formats. The missing data is filled by predefined format conversion rules and mean or interpolation methods to ensure the integrity and consistency of the data. The GeoPandas library is used to read the geological structure file and perform spatial connection to accurately match the well location data with the geological structure attributes. In addition, the well group is reasonably divided and the stratigraphic information is accurately matched to ensure that the data at each depth point contains the corresponding geological structure attributes. This innovative preprocessing method improves the data quality and lays a solid foundation for subsequent modeling.

[0020] 2. Multi-dimensional data enhancement

[0021] The present invention innovatively combines multi-dimensional data enhancement technologies such as spatial enhancement, feature enhancement, and depth enhancement. Virtual well data that conforms to geological laws is generated through geological simulation and Kriging interpolation, which improves the model's adaptability to sparse data and complex geological conditions. At the same time, adding Gaussian noise and smoothing methods enhances the robustness of the data and ensures the stability and reliability of the model in practical applications. These data enhancement methods not only expand the diversity of the data set, but also ensure the geological rationality and effectiveness of the generated data.

[0022] 3. Dynamic adjacency matrix construction

[0023] The present invention proposes a method for constructing a dynamic adjacency matrix based on the geographical location of the wells, geological similarity, and depth sequence similarity. By calculating the spatial distance, lithology similarity score, and dynamic time warping (DTW) distance, the spatial and geological relationships between wells are accurately captured, and a spatial adjacency matrix, a geological similarity adjacency matrix, and a depth sequence adjacency matrix are generated. This method effectively optimizes the effect of feature propagation and fusion in graph convolutional networks, and improves the model's ability to understand and process complex geological data.

[0024] 4. Adaptive multimodal feature fusion model

[0025] The present invention designs an adaptive multimodal feature fusion model, which comprehensively uses graph convolutional networks (GCN), multi-head attention mechanisms and gating mechanisms to achieve efficient integration of multiple adjacency matrix information. By sharing the feature extraction layer, the model can uniformly process well logging data and geological structure information, extract high-level feature representations containing temporal and spatial correlations, avoid repeated feature learning, and reduce computational complexity. At the same time, the multi-head attention mechanism and gating mechanism enhance the selectivity and expressiveness of features, further improving the prediction performance and computational efficiency of the model.

[0026] 5. Multi-task joint learning and physical constraint regularization

[0027] The present invention introduces a multi-task joint learning framework, which realizes synchronous high-precision prediction of PEM, POR, SOI and lithology classification by sharing feature extraction layer, multi-task loss balance module and multi-task prediction module. The adaptive loss weight adjustment mechanism dynamically optimizes the weight distribution according to the loss changes of each task, ensures the balanced optimization of different tasks, and prevents a certain task from dominating the overall optimization process, thereby comprehensively improving the prediction accuracy of all tasks. At the same time, the physical constraint regularization module incorporates geophysical laws into model training, and ensures the geological rationality and reliability of the prediction results by constraining the relationships between porosity and depth, water saturation and depth, porosity and density, and water saturation and relative permeability. These innovative designs not only improve the comprehensive performance of multi-task prediction, but also ensure the scientificity and reliability of the prediction results in practical applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 is a flow chart of the method of the present invention;

[0029] Figure 2 This is the structural diagram of the main parts of the model. DETAILED DESCRIPTION

[0030] The technical solutions in the embodiments of the present invention are described clearly and completely below. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0031] Embodiment 1

[0032] A well logging interpretation method based on adaptive multi-modal feature fusion of multi-task joint learning specifically includes the following steps:

[0033] S1: Data input and preprocessing stage: It is necessary to obtain logging data of multiple well groups in different blocks, and the logging data include Depth (depth), AC (acoustic time difference), CAL (well diameter), CN (neutron), CON1 (induction logging), DEN (density), EPOR (effective porosity), GR (gamma ray), K (potassium content), KTH (uranium-free gamma), PERM (absolute permeability), RS (shallow lateral resistivity), NMR (nuclear magnetic resonance logging), SH (shale content), SP (spontaneous potential), TH (thorium content), U (uranium content), UPOR (total porosity), interpretation results (permeability PEM, porosity POR, formation saturation SOI), well location data (including block ID, well ID, logging depth and latitude and longitude coordinates of the well), geological structure data (top surface depth, bottom surface depth, formation thickness, structural elevation, formation number, lithology code). The lithology code is a four-digit code, the first two digits represent the main lithology and the last two digits represent the secondary lithology.

[0034] In this process, the common problems encountered are inconsistent data formats or missing data. In order to solve these problems, the present invention adopts predefined format conversion rules and uses the mean or interpolation method to fill the missing data, thereby ensuring the integrity and consistency of the input data.

[0035] Use the geopandas library in Python to read the geological structure file and store it in GeoDataFrame. Since the well location information does not contain coordinate system information, it is necessary to create geometric points using the longitude and latitude of the well location data, and then load the well location points into the geological structure model (GeoDataFrame) and specify the coordinate system.

[0036] If the coordinate systems of the well location data and geological structure data are different, coordinate conversion is required.

[0037] The well location and geological structure data are spatially joined to obtain the geological structure attributes at the well location. If the well location point intersects with the geological structure polygon (for points and polygons, the point is actually located inside the polygon), the corresponding geological structure attributes are attached to the well location data.

[0038] The joined well location and geological structure attribute association data contains all columns from the well location data and the geological structure data.

[0039] The data associated with well locations and geological structure attributes and logging data are first classified according to block ID, and then different well groups in the same area are preliminarily divided according to the geological structure map taking into account the influence of faults, folds, etc., and then merged and saved according to the well ID of each well group. After the merger, the data of each depth point contains the geological structure attributes (stratigraphic number, lithology) of the corresponding well.

[0040] For each depth point in each well, determine whether its depth is within the range of a specific formation (based on the top and bottom depths of the formation). If a well has multiple formation segments, you need to ensure that the formation information is correctly matched. If each well has only one target formation, it can be used directly. If each well has multiple formation segments, you need to include information about each formation segment in the geological structure attributes.

[0041] For each depth point, according to the depth of the depth point, the stratigraphic section of the well is traversed to find the matching stratigraphic unit and return the corresponding geological structure attributes. After finding the stratigraphic unit to which each depth point belongs, the geological structure attributes of the stratigraphic unit are attached to the data of the depth point. In this way, the splicing of well logging feature data and geological structure data is realized, which is convenient for subsequent use.

[0042] Finally, the eigenvalues ​​and labels are divided. The eigenvalues ​​include the logging data of each well: Depth, AC, CAL, CN, CON1, DEN, EPOR, GR, K, KTH, PERM, RS, NMR, SH, SP, TH, U, UPOR, the longitude coordinates of the well, the latitude coordinates of the well, and the structural elevation, a total of 23 dimensions. The labels are: well permeability, porosity, formation saturation, and lithology, a total of 4 dimensions.

[0043] S2: Perform hybrid enhancement on the input data: Use hybrid data augmentation methods to diversify the data in spatial, feature, and depth dimensions to improve the robustness of the model.

[0044] In the data enhancement step, the present invention adopts a multi-dimensional hybrid data enhancement method, including space enhancement, feature enhancement and depth enhancement, to address the problem of uneven distribution of well data and improve the model's adaptability to sparse data.

[0045] The spatial enhancement method uses geological modeling and simulation technology to simulate the adjacency relationship between different wells and generate virtual data samples with various spatial structures that conform to geological laws. During the operation, a geological model is constructed based on the geological structure and logging characteristics of existing wells to simulate the spatial distribution and attribute changes of the strata. On this basis, new virtual well data are generated. These virtual wells maintain a reasonable spatial relationship with the original wells in terms of geographical location, and their geological characteristics are consistent with the surrounding strata.

[0046] First, based on the merged data of S1, the geological structure, stratigraphic sequence and logging characteristics of the original well are analyzed to determine the key geological parameters and change trends. Then, the geological simulation software is used to determine the spatial scope of the geological model according to the study area and the location of the existing wells, and a stratigraphic framework model is established based on the stratigraphic sequence to define the spatial distribution and extension of each layer. Then, the lithological parameters (such as porosity, permeability, etc.) in the logging data are assigned to the grid of the geological model, and geological structural elements such as faults and folds are incorporated into the model to ensure that the model reflects the real geological structural characteristics. Appropriate locations are selected in the geological model to generate data for virtual wells. These locations should meet the following requirements:

[0047] (1) Geological continuity: The location of the virtual well should be within the effective range of the geological model and have geological continuity with the original well, avoiding crossing faults or discontinuous geological interfaces.

[0048] (2) Spatial rationality: The distance between the virtual well and the original well should be within a reasonable range, usually set within 500 meters, to ensure spatial correlation.

[0049] (3) Similarity of geological characteristics: Select a location with similar geological characteristics to the original well to ensure that the data of the virtual well is consistent with the actual geological conditions.

[0050] (Method 2) Based on the existing well location data, generate virtual well locations that meet geological continuity and spatial rationality. Construct a KD tree based on the longitude and latitude coordinates of the existing well groups in each stratum of each block to facilitate the search for the nearest neighbor wells. The number of virtual wells generated by the same well group is set to 0-10. The longitude and latitude range for generating virtual wells is obtained by calculating the maximum and minimum values ​​of the existing well coordinates in the same stratum. A longitude and latitude is randomly generated by looping to generate the required number of virtual wells. The generated virtual well location is then queried through the KD tree to find the original well closest to the location and calculate the Euclidean distance between the two. If the distance is less than or equal to 0.02 degrees (approximately equal to 2000m) and greater than or equal to 0.003 degrees (approximately equal to 300m), it is considered that the location of the virtual well is within a reasonable spatial range and meets the requirements of geological continuity. Therefore, the virtual well location data and virtual well ID are stored in the list of the stratum of the block.

[0051] To ensure the reliability of virtual well data, the generated logging features should be based on the data of surrounding wells. The present invention adopts ordinary Kriging interpolation method, uses the logging data and spatial position of existing wells, and estimates the geological attribute values ​​at the unsampled position (virtual well) according to the variation function model.

[0052] Prepare the grid coordinates required for interpolation according to the latitude and longitude data of the virtual well, and use gridx and gridy to extract the longitude and latitude of the virtual well as the interpolation target point. When generating the 18-dimensional eigenvalues ​​of the logging data of the virtual well, first set a Kriging interpolation object and pass in the eigenvalue of the original well and the corresponding longitude and latitude coordinates lon and lat, and specify the variogram model as linear, interpolate the position of the virtual well, generate the predicted value and variance, and finally obtain the 18-dimensional eigenvalues ​​of the logging data of the virtual well.

[0053] For lithology data, since it is a categorical variable, Kriging interpolation cannot be used directly. The nearest neighbor method is used for lithology prediction. By constructing a KD tree, the 1-3 nearest original wells near the location of each virtual well are queried to obtain the lithology information and formation information of these wells. Then, based on the frequency of occurrence of each formation lithology category in these 1-3 nearest wells, the lithology category with the most occurrences in the same formation is copied to the virtual well. The virtual data generated in this way not only retains spatial rationality, but also conforms to the law of change of geological characteristics.

[0054] Through this spatial enhancement method based on geological simulation, the spatial diversity of the data set is expanded, and the adaptability of the model to different spatial structures and geological conditions is improved. At the same time, the geological irrationality that may be caused by directly adjusting the well location is avoided, ensuring the scientificity and effectiveness of data enhancement.

[0055] The feature enhancement method adds a small amount of random noise to the logging features, making the generated samples slightly different from the original data, thereby improving the model's ability to resist noise interference. According to the changing rules of different eigenvalues, a Gaussian noise with a mean of 0 and a variance of a certain proportion of its standard deviation (such as 1% or 5%) is added to the features of each logging data, so that the added noise amplitude matches the numerical scale of the feature to generate reasonable virtual data samples, thereby avoiding the introduction of unreasonable outliers due to enhancement.

[0056] For depth enhancement, the present invention adopts nonlinear interpolation method and sliding window technology to generate more depth segment samples and improve the model's adaptability to depth dimension data. In the interpolation of the depth sequence of each well, the spline interpolation method is used to more accurately fit the nonlinear changes of geological characteristics with depth. Spline interpolation constructs a smooth, polynomial function between known data points, making the interpolation curve continuous at each data point and having a continuous derivative, thereby better simulating the changing trend of geological attributes.

[0057] During the operation, the logging characteristic data of adjacent depth points are first collected. i and d i+1 The remaining 17-dimensional logging eigenvalues ​​f(d i ) and f(d i+1 ). Using the known logging eigenvalues, construct the spline function S(d) and i With d i+1 Generate additional interpolation points d between j , calculate the corresponding logging parameter value f(d j )=S(d j ). In this way, the missing data in depth are supplemented, making the logging data more continuous and smooth in the depth dimension.

[0058] In the application of sliding window technology, the window size is set to 10 frames, that is, each window contains 10 consecutive depth point data. The sliding window moves across the entire depth range with a step size of 3 frames, extracting data fragments in each window to form a new sample. This method can capture the local changes of geological features at depth, extract deep dynamic features, and enhance the model's ability to learn the laws of temporal changes.

[0059] To ensure the rationality and effectiveness of the generated enhanced data, it is necessary to refer to the existing geological empirical formulas or stratigraphic division results to guide the interpolation process during the interpolation and windowing process. In areas where stratigraphic interfaces or geological properties change dramatically, the interpolation method or density should be appropriately adjusted to avoid introducing unreasonable samples.

[0060] The generated enhanced data needs to be quality verified. Through statistical analysis, the distribution characteristics of the enhanced data and the original data are compared to ensure that they are physically reasonable and do not introduce outliers or data that violate geological laws. If unreasonable data points are found, they should be corrected or eliminated.

[0061] Finally, there is noise processing. Since there may be acquisition errors and environmental noise in actual logging data, these noises will affect the prediction accuracy of the model. In the noise processing part, the present invention adopts a variety of methods to detect and process acquisition errors and environmental noise in logging data to ensure data quality. Specific methods include statistical outlier detection, mean interpolation substitution, high-frequency noise smoothing, and robust smoothing filtering.

[0062] A statistical outlier detection method is used to identify outliers in the data that are significantly deviated from the normal range. Specifically, for each logging feature except depth, the mean (μ) and standard deviation (σ) of the feature are calculated, and a detection threshold is set to identify outliers. The threshold is set to 3 times the standard deviation. If a data point x i Deviating from the mean by more than 3 standard deviations (|x i -μ|>3·σ), then mark the point as an outlier. This method can effectively detect statistically significant abnormal data.

[0063] After detecting the outlier, traverse the 17-dimensional logging feature value data again to obtain the outlier mask and the row index, and determine whether the outlier is a boundary point. If the outlier is not a boundary point, replace the point with the average value of its adjacent points before and after. The formula is:

[0064] If the point is a boundary point, that is, the data is located in the first row or the last row, the mean of the feature is directly used to replace the outlier.

[0065] In addition, in order to reduce high-frequency noise in the feature data, the present invention also introduces smoothing processing in the data processing stage. The sliding average method divides the data of the logging features other than PERM into segments and averages the data points in each segment to smooth out subtle fluctuations. The specific formula is:

[0066]

[0067] Where n is the window size, and the width of the sliding window is adjusted according to the noise level. When the window is large, the impact of high-frequency noise can be reduced more smoothly.

[0068] Traverse the logging eigenvalues, set the sliding window size to 5 data points, the window is at the center of the current point, and calculate the average value of the data points in the window as the processed eigenvalue data.

[0069] For data segments with particularly obvious noise (PERM data segments), the present invention uses the median filtering method to replace the value of each data point with the median of the data points in its neighborhood, avoiding the influence of outliers on the smoothing effect. The specific formula is:

[0070]

[0071] The window size is k, and the median of all data points in the window is used as the value of the position.

[0072] Set its window size to 5 data points, with the window centered at the current point.

[0073] Next is data normalization. Since the data magnitudes of different logging characteristics vary greatly, if they are directly input into the model, the weights of certain features may be too large or too small, thereby affecting the stability and convergence speed of the training. Therefore, the present invention adopts a normalization method to convert the data ranges of all features into the same scale, thereby ensuring that the feature value size does not affect the learning effect of the model. Specifically, two methods, Min-Max normalization and Z-score standardization, are used. For data without obvious outliers such as (AC, UPOR, SH, CAL, CN, DEN, etc.), Min-Max normalization is used to map the data to the [0,1] interval, and the calculation formula is:

[0074]

[0075] This method can ensure that the distribution range of the data is within a fixed interval, which is convenient for model processing. In addition, for data with extreme values ​​(PERM, GR), Z-score standardization is used to further reduce the impact of extreme values ​​by adjusting the data to a distribution with a mean of 0 and a variance of 1.

[0076] S3: Dynamic adjacency matrix construction module: Generate spatial adjacency matrix through the geographic spatial location of wells; generate geological similarity adjacency matrix based on the similarity of geological characteristics between wells; apply sliding depth window to generate depth sequence adjacency matrix.

[0077] 301. Generate a spatial adjacency matrix. In order to capture the geographic adjacency relationship between wells, the two-dimensional spatial distance between wells is calculated based on the well location information (well ID, longitude, latitude) of the same block. The specific calculation formula is:

[0078]

[0079] Among them, x, y, and z are the longitude, latitude, and depth of the well location, respectively. Then, the distance threshold is set to 0.005 degrees (about 500 meters), and only the well pairs with a distance less than or equal to the threshold are retained as adjacency relationships. For the well pairs that meet the distance requirements, self-connection is avoided and the distance is inversely proportional (such as ) Calculate the edge weights and generate a spatial adjacency matrix. The weights exceeding the threshold are set to 0. At the same time, if the distance of the same well is 0, the weight will be infinite, so the weight is also 0 when the distance is set to 0.

[0080] Lithology similarity is introduced into the spatial adjacency matrix. For each pair of wells in the same block, according to the lithology coding, the first two digits represent the main lithology and the last two digits represent the secondary lithology.

[0081] In order to speed up the calculation, the lithology similarity score can be calculated in parallel.

[0082] By establishing a process pool, an empty list pool_args is initialized to store all parameter pairs that need to be processed in parallel. Through nested double loops, all formation lithology pairs are traversed (ensuring that each pair of wells is calculated only once and maintaining symmetry, i.e., i <= j)), the well ID and lithology data of each pair of wells are extracted and added to pool_args, and then tasks are assigned to each process through the pool_args list to parallelly calculate the similarity scores of all formation lithologies, where the similarity scores are as follows: (1) The same primary lithology and the same secondary lithology: 1.0

[0083] (2) The primary lithology is the same, and the secondary lithology is different but similar: 0.7

[0084] (3) The primary lithology is the same, but the secondary lithology is different and dissimilar: 0.5

[0085] (4) Different primary lithology and same secondary lithology: 0.5

[0086] (5) The primary lithology is different, and the secondary lithology is different but similar: 0.3

[0087] (6) The primary lithology is different, and the secondary lithology is different and dissimilar: 0.0

[0088] Then, the lithology of multiple formations is weighted averaged according to the formation thickness to obtain the final formation lithology score of the well pair, and finally the calculated weighted average score is assigned to the corresponding well pair in the similarity matrix.

[0089] Then, based on the adjacency matrix obtained by considering the distance before, the obtained lithology similarity matrix is ​​multiplied with the previous spatial adjacency matrix to update the spatial adjacency matrix.

[0090] 302. Construct a geological similarity adjacency matrix. This matrix determines the adjacency relationship by the similarity of the 17-dimensional feature values ​​of the logging data except the depth. For each feature, the cosine similarity between wells is calculated:

[0091] Among them, T ik and T jk are the values ​​of well i and well j on feature k, respectively. All feature similarities are weighted and summed to obtain the comprehensive similarity. The similarity threshold is set to 0.65, well pairs with similarity higher than the threshold are retained, geological similarity adjacency relationships are established, and the similarity values ​​are filled into the adjacency matrix as edge weights.

[0092] 303. Generate a depth sequence adjacency matrix. In order to capture the dynamic correlation of well logging data, the gradually increasing depth in the logging data is regarded as pseudo time, and a sliding depth window is used to process the depth sequence data.

[0093] The well logging feature data divided and stored according to S2 (the depth window size is set to 10 and the sliding step size is 3)

[0094] The dynamic depth warping (DTW) distance of the well pair data is calculated in each depth window in a multi-process parallel manner to measure the similarity of the depth sequence. The specific formula is:

[0095]

[0096] The DTW distance calculation function is used to calculate the DTW distance of all sliding window pairs between each well pair. The smaller the distance, the more similar the two windows are. The distance of each pair of windows is accumulated. If a window pair exists, its average distance is returned; otherwise, infinity is returned (indicating that it cannot be compared). Finally, the average DTW distance of all window pairs between the two wells is obtained.

[0097] The threshold is dynamically adjusted according to the distribution of DTW distance, the DTW distance threshold is set to 95% of the distance, and then the adjacency matrix is ​​constructed according to the threshold.

[0098] Finally, the three adjacency matrices of each block are calculated through a loop, and then the corresponding adjacency matrices of each block are spliced ​​on the diagonal, and the adjacency relationship between blocks is set to 0, and finally three global adjacency matrices are obtained.

[0099] S4: Adaptive multimodal feature fusion model construction:

[0100] (1) Shared feature extraction layer

[0101] The well logging features and geological structure features are spliced ​​in the feature dimension to form the initial feature matrix X. The well logging features of each well have 17 dimensions and the geological structure features have 6 dimensions. The dimension of the spliced ​​input feature matrix X is N×23 (N is the total number of wells in all blocks).

[0102] X is transformed using a linear layer to map the feature dimension from 23 to a hidden layer dimension of 128. The transformed features are activated by the ReLU nonlinear activation function to obtain the initial feature representation H0, whose dimension is N×128.

[0103] The three adjacency matrices obtained by S3 are spatial adjacency matrices A. space , geological similarity adjacency matrix A geological and the deep sequence adjacency matrix A depth, the dimension of each adjacency matrix is ​​N×N. Give each adjacency matrix a learnable weight coefficient α i (i=1,2,3), the initial value is set to Ensure that each adjacency matrix contributes equally to feature extraction in the initial stage. Normalize each adjacency matrix to obtain

[0104] Array The calculation formula is as follows:

[0105]

[0106] Among them, D i is the adjacency matrix A i The degree matrix of .

[0107] Use the Softmax function to normalize the weight coefficient and get the normalized weight β i :

[0108]

[0109] Comprehensive adjacency matrix It is obtained by weighted summation of three normalized adjacency matrices:

[0110] By weighted summing of three normalized adjacency matrices, we get:

[0111]

[0112] Using the comprehensive adjacency matrix Perform graph convolution on the initial feature representation H0. Multiply to get the propagated features Its dimension is N×128.

[0113] Then a linear transformation layer is used to map the feature dimension from 128 to 64, and a ReLU activation function is applied to generate the graph convolutional layer output H1, whose dimension is N×64.

[0114] In order to enhance the selectivity of the model for important features, a gating mechanism is introduced in step S4. The output H1 of the graph convolutional layer passes through a linear layer and a Sigmoid activation function to generate a gating coefficient matrix G, whose dimension is also N×64 and the value range is between 0 and 1. By performing an element-wise product operation on H1 and g, the gated feature representation H is obtained. gated , with a dimension of N × 64. This process enables the model to dynamically adjust the transfer weights of each feature channel, highlight key features, and suppress redundant information.

[0115] After the gating mechanism, the model further integrates an adaptive adjacency matrix attention mechanism to optimize the weight distribution of different adjacency matrices. The model defines a new learnable parameter γ for each adjacency matrix i (i=1,2,3), the initial value is set to Through the Softmax function, γ i Normalized to weight δ i , ensuring that the weight sum is 1. Comprehensive adjacency matrix It is calculated by the following formula:

[0116]

[0117] Among them, A i is the original adjacency matrix.

[0118] Then the multi-head attention mechanism is introduced to further improve the expressiveness of features.

[0119] The number of attention heads is set to 4, and each attention head processes 16 subspaces of the feature dimension. For each attention head, the model first generates a query vector, a key vector, and a value vector through linear transformation, and then calculates the attention coefficient between nodes. Based on these coefficients, the features of neighboring nodes are weighted and summed to obtain the new feature representation of each node under this attention head. The outputs of all attention heads are concatenated in the feature dimension to form the final attention output H attn , whose dimension is N×64.

[0120] In order to integrate the gated features H gated And the output H of the multi-head attention mechanism attn The model uses a residual connection to add the two and pass the ReLU activation function to generate the final fusion feature representation H final , whose dimension is N×64.

[0121] Finally, for different prediction targets, the model designs independent multi-task output layers. For lithology prediction (classification task), a fully connected layer is used to combine the fusion feature H final Mapped to the number of categories C, and applying the Softmax activation function, the probability distribution y of the lithology category is output lithilogy , whose dimension is N×C. For the permeability y perm , porosity y por and formation saturation y soi prediction (regression task), using independent fully connected layers to fuse the fusion features H final Mapped to a single scalar value, the output dimensions are all N×1.

[0122] S5 multi-task joint learning: mainly through the shared feature extraction layer, multi-task loss balance module and multi-task prediction module to achieve the joint prediction of three well logging indicators and lithology classification of PEM, POR and SOI. The goal of this module is to effectively utilize the shared spatial features and improve the comprehensive prediction accuracy of multiple tasks.

[0123] 501. In the shared feature extraction layer, the model sets a shared feature extraction layer for different tasks. The shared feature extraction layer is transformed from H final A higher-level feature representation is extracted, which not only contains the basic characteristics of well logging indicators such as PEM, POR, SOI, but also incorporates the geological information required for lithology classification. By sharing the feature extraction layer, the model avoids learning feature representations for each task separately, reducing the computational burden, while ensuring the consistency and complementarity of features between tasks.

[0124] 502. Use adaptive adjustment of the loss weight of each task to ensure balanced optimization of each task during model training. Since the loss scale and learning rate of each task are different, directly adding their loss functions may cause the loss of one task to dominate the overall optimization process, resulting in a decrease in the prediction accuracy of the remaining tasks. To solve this problem, a weight factor α can be introduced, and the loss function L of each task task is adaptively adjusted, the formula is:

[0125] L total =α1L PEM +α2L POR +α3L SOI +α4L lithology

[0126] Among them, L PEM , L POR , L SOI are the mean square error losses of permeability, porosity and formation saturation prediction, L lithology is the cross entropy loss for lithology classification. The weight factors α1, α2, and α3 are dynamically adjusted during the model training process to achieve adaptive weight allocation for each task.

[0127] Adaptive weight allocation uses a dynamic weighting method to enable the model to find the best weight configuration during the continuous optimization process, ensuring the reasonable contribution of each task's loss to the total loss and avoiding excessive impact of a certain task on the overall training. By introducing an automatic weight adjustment mechanism, αl, α2, and α3 are dynamically adjusted according to the loss changes of each task, or the weights are adjusted based on the uncertainty of the task, so that the model can achieve balanced optimization in the multi-task learning process.

[0128] Finally, the model outputs to the multi-task prediction module, which contains independent output layers for multiple tasks to generate the final prediction values ​​and lithology classification results of PEM, POR and S0I. For each task, the output layer contains a fully connected layer to map the shared features to the target logging indicators. The output layer of each task is respectively passed through the corresponding loss function L PEM , L POR , L SOI and L lithology Make optimizations.

[0129] For each task, an independent fully connected layer is designed to share the shared features H shared Map to the target logging indicator or category. The specific implementation is as follows:

[0130] Permeability (PEM) prediction: A fully connected layer maps the shared features to a single scalar value and outputs the predicted value of the permeability

[0131]

[0132] Porosity (POR) prediction: A fully connected layer maps the shared features to a single scalar value and outputs the predicted porosity value

[0133]

[0134] Formation saturation (SOI) prediction: A fully connected layer maps the shared features to a single scalar value and outputs the predicted value of formation saturation

[0135]

[0136] Lithology classification: A fully connected layer maps the shared features to the number of categories C, and applies the Softmax activation function to output the probability distribution of lithology categories. where CCC is the number of lithology categories.

[0137]

[0138] Through the multi-task prediction module, the model can independently perform regression prediction and lithology classification of PEM, POR, and SOI based on shared features, thereby achieving more accurate multi-task logging index prediction. The independent output layer of each task ensures that each task performs personalized prediction based on shared features, avoiding interference between tasks and improving the overall prediction performance.

[0139] Training of S6 model:

[0140] Using the preprocessed and enhanced feature matrix X (dimension is N×17) obtained in steps S1 and S2, the adjacency matrix set [A space ,A geological ,A depth ] (the dimension of each adjacency matrix is ​​N×N), and the corresponding label data y lithology ,y perm ,y porosity ,y saturation ,The data set is evenly sampled in each block according to the blocks, and the data set is divided into training set, validation set, and test set, with a ratio of 8:1:1.

[0141] Then define the loss function of the classification task (lithology prediction) as cross entropy loss (CrossEntropyLoss),

[0142] The loss function for the regression tasks (PEM, POR, SOI) is defined as the mean square error loss (MSELoss).

[0143] The model built using S4 is trained in a sliding window loop with a size of 10 and a step size of 3, and the cumulative loss is initialized to 0.

[0144] Then the feature data, three adjacency matrices, and four label data of the N wells processed each time are moved to the designated computing device to enable batch processing.

[0145] Then clear the accumulated gradients in the optimizer to avoid the impact of gradient accumulation on the current step.

[0146] The input feature X and the three adjacency matrices are input into the model to obtain the prediction output of each task, including the prediction results of lithology classification (N, C), permeability (N, 1), and porosity (N, 1).

[0147] Then the lithology predicted losses and the predicted losses of PEM, POR, and SOI are calculated.

[0148] Then the physical constraint regularization loss is calculated to ensure that the model predictions conform to actual geological laws:

[0149] 601. Empirical Relationship between Porosity and Depth

[0150] In reservoir geology, porosity usually decreases with increasing depth, a phenomenon mainly due to the formation compaction effect. A common empirical formula used to describe the relationship between porosity and depth is as follows:

[0151] φ=φ0e -c·d

[0152] Where: φ represents porosity, φ0 represents the initial porosity at the shallow layer, c is the compaction coefficient (depending on the formation properties, usually between 0.0001 and 0.0005), and d represents the depth.

[0153] This formula can be used as a regularization constraint to make the relationship between the porosity and depth predicted by the model consistent with the compaction effect. The following regularization term is added to the loss function to calculate the mean square error between the predicted porosity and the expected porosity:

[0154]

[0155] 602. Relationship between water saturation and depth:

[0156] Water saturation S w It usually decreases with increasing depth, especially in deeper parts of the reservoir, where the water saturation tends to be lower due to capillary pressure. According to the research of JD Buckley et al. (1995), the following empirical formula can be used to describe the relationship between depth and water saturation:

[0157] S w =S wc +(S wi -S wc ) -α·d

[0158] Where: S w c is the residual water saturation, which represents the final water saturation value at a deeper level, S wi is the shallow water saturation value, α is the water saturation variation coefficient (generally between 0.0001 and 0.0003), and d represents the depth.

[0159] This formula can be used as a regularization constraint to ensure that the depth variation trend of water saturation is maintained. The regularization term is as follows, which is used to calculate the mean square error between the predicted water saturation and the expected water saturation:

[0160]

[0161] 603. Relationship between porosity and density:

[0162] Formation density ρ usually increases with depth, while porosity decreases with increasing density. The commonly used relationship is as follows

[0163]

[0164] Where: matrix is the matrix density of the formation rock, which is usually a known value (for example, the matrix density of limestone is about 2.71 g / cm 3 );ρ fluidis the density of the fluid (such as water or oil, about 1g / cm 3 ).

[0165] This formula can express the relationship between porosity and density as a regularization constraint:

[0166]

[0167] 604. Relationship between water saturation and relative permeability

[0168] In reservoir description, water saturation S w and relative permeability k r There is a correlation between

[0169] k r =k r0 ·(1-S w ) n

[0170] Where: k r0 is the initial permeability when the saturation is zero, and n is the permeability curve index, which usually takes a value of 2 to 4, depending on the reservoir type.

[0171] This relationship can be used as a constraint in the model to ensure reasonable water saturation and permeability:

[0172]

[0173] 605. Complete loss function:

[0174] In order to implement these regularization constraints, the above regularization terms can be added to the calculation of the loss function and the model parameters can be updated through back propagation. The overall form of the loss function is:

[0175] L=L base +λ1R porosity-depth +λ2R water-saturation-depth +λ3R porosity-density

[0176] +λ4R permeability-saturation

[0177] L base =L lithology +L perm +L porosity +L saturation

[0178] The gradient of the total loss L with respect to the model parameters is then calculated through backpropagation.

[0179] The Adam optimizer is used to update the model parameters according to the calculated gradients to minimize the total loss.

[0180] Add the total loss of the current batch to the total loss, and then calculate the average loss of the current epoch

[0181] After each epoch, the model is evaluated using the validation set to calculate the performance indicators RMSE and accuracy of each task.

[0182] Finally, according to the performance of the validation set, the loss weight α is dynamically adjusted i and the regularization weight β i .

[0183] Embodiment 2

[0184] The model uses a combination of graph convolutional networks (GCN), multi-head attention, gating mechanism and adaptive loss weight adjustment algorithm to achieve high-precision joint prediction of permeability (REM), porosity (POR), formation saturation (SOI) and lithology classification. By dynamically constructing and integrating spatial adjacency matrix, geological similarity adjacency matrix and depth sequence adjacency matrix, the model can effectively capture the complex spatial and geological relationships between wells. The physical constraint regularization module integrates geophysical laws into the model training process to ensure that the prediction results are consistent with the actual geophysical laws on the basis of statistical accuracy. Through advanced machine learning algorithms and specific graph neural network structure design, this technical solution significantly improves the comprehensive prediction ability of multi-well logging data and the robustness of the model. It is widely used in oil and gas reservoir evaluation, geological modeling and related geological data analysis, and greatly promotes the effective development and management of oil and gas resources.

[0185] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A well logging interpretation method based on adaptive multimodal feature fusion based on multi-task joint learning, characterized in that: The following steps are involved: S1: Data input and preprocessing, including obtaining multi-well group logging data, well location data and geological structure data, and performing format conversion, missing value filling, coordinate conversion and spatial connection on the data to ensure data integrity and consistency; S2: Data enhancement processing, using a multi-dimensional hybrid data enhancement method, including spatial enhancement, feature enhancement, and depth enhancement, to generate virtual well data that conforms to geological laws, improving the model's robustness and adaptability to sparse data; S3: Dynamic adjacency matrix construction, based on the geographic spatial location, geological characteristics similarity and depth sequence similarity of the wells, the spatial adjacency matrix, geological similarity adjacency matrix and depth sequence adjacency matrix are constructed respectively to accurately capture the complex spatial and geological relationships between wells; S4: Adaptive multimodal feature fusion model construction, using graph convolution, multi-head attention mechanism and gating mechanism to extract high-level features from the fused features, and integrate the information of multiple adjacency matrices through dynamic weighting to generate a unified feature representation; S106 multi-task joint learning, through the shared feature extraction layer, multi-task loss balance module and multi-task prediction module, achieves simultaneous high-precision prediction of PEM, POR, SOI and lithology classification; S5: Model training, using the preprocessed and enhanced feature matrix and the constructed adjacency matrix to train the constructed multimodal feature fusion model, adopting a multi-task joint learning framework, combining cross entropy loss and mean square error loss with physical constraint regularization terms to ensure that the prediction results conform to the laws of geophysics.

2. The well logging interpretation method based on multi-task joint learning and adaptive multi-modal feature fusion according to claim 1, characterized in that: In S1: Data input and preprocessing stage: It is necessary to obtain logging data of multiple well groups in different blocks, including depth, acoustic time difference, wellbore, neutron, induction logging, density, effective porosity, gamma ray, potassium content, uranium-free gamma, absolute permeability, shallow lateral resistivity, nuclear magnetic resonance logging, shale content, natural potential, thorium content, uranium content, total porosity, interpretation results (permeability PEM, porosity POR, formation saturation SOI), well location data (including block ID, well ID, logging depth and latitude and longitude coordinates of the well), geological structure data (top depth, bottom depth, formation thickness, structural elevation, formation number, lithology code); Predefined format conversion rules are used and missing data are filled using mean or interpolation methods to ensure the integrity and consistency of input data; Use the geopandas library in Python to read the geological structure file and store it in GeoDataFrame; since the well location information does not contain coordinate system information, it is necessary to create geometric points using the longitude and latitude of the well location data, and then load the well location points into the geological structure model and specify the coordinate system; If the coordinate systems of the well location data and geological structure data are different, coordinate conversion is required; The well location is spatially connected with the geological structure data to obtain the geological structure attributes at the well location. If the well location intersects with the geological structure polygon, the corresponding geological structure attributes are added to the well location data. The connected well location and geological structure attribute association data contains all columns from the well location data and the geological structure data; The data associated with well locations and geological structure attributes and logging data are first classified according to the block ID, and then different well groups in the same area are preliminarily divided according to the geological structure map, taking into account the influence of faults, folds, etc., and then merged and saved according to the well ID of each well group. After merging, the data of each depth point contains the geological structure attributes (stratum number, lithology) of the corresponding well; For each depth point of each well, determine whether its depth is within the range of a specific formation (based on the top and bottom depths of the formation); if a well has multiple formation segments, it is necessary to ensure that the formation information is correctly matched; if each well has only one target formation, it can be used directly; if each well has multiple formation segments, it is necessary to include information about each formation segment in the geological structure attributes; For each depth point, according to the depth of the depth point, the stratigraphic section of the well is traversed to find the matching stratigraphic unit and return the corresponding geological structure attributes; after finding the stratigraphic unit to which each depth point belongs, the geological structure attributes of the stratigraphic unit are attached to the data of the depth point; in this way, the splicing of the well logging feature data and the geological structure data is realized, which is convenient for subsequent use; Finally, the eigenvalues ​​and labels are divided. The eigenvalues ​​include the logging data of each well: Depth, AC, CAL, CN, CON1, DEN, EPOR, GR, K, KTH, PERM, RS, NMR, SH, SP, TH, U, UPOR, the longitude coordinates of the well, the latitude coordinates of the well, and the structural elevation, a total of 23 dimensions. The labels are: well permeability, porosity, formation saturation, and lithology, a total of 4 dimensions.

3. The well logging interpretation method based on multi-task joint learning and adaptive multi-modal feature fusion according to claim 1, characterized in that: In step S2: hybrid enhancement processing is performed on the input data: a hybrid data enhancement method is used to diversify the data in space, feature, and depth dimensions to improve the robustness of the model; In the data enhancement step, the present invention adopts a multi-dimensional hybrid data enhancement method, including space enhancement, feature enhancement and depth enhancement, to address the problem of uneven distribution of well data and improve the model's adaptability to sparse data.

4. The well logging interpretation method based on adaptive multi-modal feature fusion based on multi-task joint learning according to claim 3 is characterized by: The spatial enhancement method is: using geological modeling and simulation technology to simulate the adjacency relationship between different wells, and generate virtual data samples with various spatial structures that conform to geological laws; during the operation, based on the geological structure and logging characteristics of existing wells, a geological model is constructed to simulate the spatial distribution and attribute changes of the strata; on this basis, new virtual well data is generated. These virtual wells maintain a reasonable spatial relationship with the original wells in terms of geographical location, and their geological characteristics are consistent with the surrounding strata; First, based on the merged data from S1, the geological structure, stratigraphic sequence and logging characteristics of the original well were analyzed to determine the key geological parameters and their changing trends; Then, using geological simulation software, according to the study area and the location of existing wells, determine the spatial scope of the geological model, and according to the stratigraphic sequence, establish a stratigraphic framework model to define the spatial distribution and extension of each layer. Then, assign the lithological parameters (such as porosity, permeability, etc.) in the logging data to the grid of the geological model, and incorporate geological structural elements such as faults and folds into the model to ensure that the model reflects the real geological structural characteristics; select appropriate locations in the geological model to generate virtual well data, and these locations should meet the following requirements: (1) Geological continuity: The location of the virtual well should be within the effective range of the geological model and have geological continuity with the original well, avoiding crossing faults or discontinuous geological interfaces; (2) Spatial rationality: The distance between the virtual well and the original well should be within a reasonable range, usually set within 500 meters, to ensure spatial correlation; (3) Similarity of geological characteristics: Select a location with similar geological characteristics to the original well to ensure that the data of the virtual well is consistent with the actual geological conditions.

5. The well logging interpretation method based on multi-task joint learning and adaptive multi-modal feature fusion according to claim 3 is characterized by: The spatial enhancement method is as follows: based on the existing well location data, a virtual well location that meets the geological continuity and spatial rationality is generated; a KD tree is constructed according to the longitude and latitude coordinate list of the existing well group in each stratum of each block to facilitate the search for the nearest neighbor well, and the number of virtual wells generated by the same well group is set to 0-10, and a longitude and latitude range for generating virtual wells is obtained by calculating the maximum and minimum values ​​of the existing well coordinates in the same stratum, and a longitude and latitude is randomly generated by looping to generate the required number of virtual wells; the generated virtual well location is then queried through the KD tree to find the original well closest to the location, and the Euclidean distance between the two is calculated; if the distance is less than or equal to 0.02 degrees and greater than or equal to 0.003 degrees, it is considered that the location of the virtual well is within a reasonable spatial range and meets the requirements of geological continuity, so the virtual well location data and virtual well ID are stored in the list of the stratum of the block; To ensure the reliability of virtual well data, the generated logging features should be based on the data of surrounding wells. The present invention adopts ordinary Kriging interpolation method, uses the logging data and spatial position of existing wells, and estimates the geological attribute values ​​at the unsampled position, i.e., the virtual well, according to the variogram model; Prepare the grid coordinates required for interpolation according to the latitude and longitude data of the virtual well, and use gridx and gridy to extract the longitude and latitude of the virtual well as the interpolation target point respectively; when generating the 18-dimensional eigenvalue of the logging data of the virtual well, first set a Kriging interpolation object and pass in the eigenvalue of the original well and the corresponding longitude and latitude coordinates lon and lat, and specify the variation function model as linear, interpolate the position of the virtual well, generate the predicted value and variance, and finally obtain the 18-dimensional eigenvalue of the logging data of the virtual well; For lithology data, since it is a categorical variable, Kriging interpolation cannot be used directly. The nearest neighbor method is used for lithology prediction. By constructing a KD tree, the 1-3 nearest original wells near each virtual well are queried to obtain the lithology information and formation information of these wells. Then, according to the frequency of occurrence of each formation lithology category in these 1-3 nearest wells, the lithology category with the most occurrences in the same formation is copied to the virtual well. The virtual data generated in this way not only retains spatial rationality, but also conforms to the law of change of geological characteristics. This spatial enhancement method based on geological simulation expands the spatial diversity of the data set and improves the model's adaptability to different spatial structures and geological conditions. At the same time, it avoids the geological irrationality that may be caused by directly adjusting the well location, ensuring the scientificity and effectiveness of data enhancement. The feature enhancement method adds a small amount of random noise to the logging features, making the generated samples slightly different from the original data, thereby improving the model's ability to resist noise interference; according to the change law of different eigenvalues, a Gaussian noise with a mean of 0 and a variance of a certain proportion of its standard deviation is added to the features of each logging data, so that the added noise amplitude matches the numerical scale of the feature to generate reasonable virtual data samples, thereby avoiding the introduction of unreasonable outliers due to enhancement; For depth enhancement, the present invention adopts nonlinear interpolation method and sliding window technology to generate more depth segment samples and improve the adaptability of the model to depth dimension data; in the interpolation of the depth sequence of each well, the spline interpolation method is used to more accurately fit the nonlinear changes of geological characteristics with depth; spline interpolation constructs a smooth, polynomial function between known data points, so that the interpolation curve is continuous at each data point and has a continuous derivative, thereby better simulating the change trend of geological attributes; During the operation, the logging characteristic data of adjacent depth points are first collected. i and d i+1 The remaining 17-dimensional logging eigenvalues ​​f(d i ) and f(d i+1 ) ; Using the known logging eigenvalues, construct the spline function S(d) and i With d i+1 Generate additional interpolation points d between j , calculate the corresponding logging parameter value f(d j )=S(d j ); In this way, the missing data in depth are supplemented, making the logging data more continuous and smooth in the depth dimension; In the application of sliding window technology, the window size is set to 10 frames, that is, each window contains 10 consecutive depth point data; the sliding window moves in the entire depth range with a step size of 3 frames, and data fragments are extracted in each window to form a new sample; This method can capture the local changes of geological features at depth, extract deep dynamic features, and enhance the model's ability to learn the laws of temporal changes; To ensure the rationality and effectiveness of the generated enhanced data, it is necessary to refer to the existing geological empirical formulas or stratigraphic division results to guide the interpolation process during the interpolation and windowing process; in areas where stratigraphic interfaces or geological attributes change dramatically, the interpolation method or density should be appropriately adjusted to avoid introducing unreasonable samples; The generated enhanced data needs to be quality verified: by comparing the distribution characteristics of the enhanced data with the original data through statistical analysis, it is ensured that it is physically reasonable and does not introduce abnormal values ​​or data that violates geological laws; if unreasonable data points are found, they should be corrected or eliminated; Finally, there is noise processing: since there may be acquisition errors and environmental noise in the actual logging data, these noises will affect the prediction accuracy of the model; in the noise processing part, the present invention adopts a variety of methods to detect and process the acquisition errors and environmental noise in the logging data to ensure data quality; the specific methods include statistical outlier detection, mean interpolation substitution, high-frequency noise smoothing processing and robust smoothing filtering, etc.; A statistical outlier detection method is used to identify outliers in the data that are significantly deviated from the normal range. Specifically, for each logging feature except depth, the mean (μ) and standard deviation (σ) of the feature are calculated, and a detection threshold is set to identify outliers. The threshold is set to 3 times the standard deviation. If a data point x i More than 3 standard deviations from the mean (|x i -μ|>3·σ), then mark the point as an outlier; this method can effectively detect statistically significant abnormal data; After detecting the outlier, traverse the 17-dimensional logging feature value data again to obtain the outlier mask and the row index, and determine whether the outlier is a boundary point. If the outlier is not a boundary point, replace the point with the average value of its adjacent points before and after. The formula is: If the point is a boundary point, that is, the data is located in the first row or the last row, the mean of the feature is directly used to replace the outlier; In addition, in order to reduce high-frequency noise in the feature data, the present invention also introduces smoothing processing in the data processing stage; the sliding average method divides the data of the logging features except PERM into segments, and averages the data points in each segment to smooth out subtle fluctuations; the specific formula is: Where n is the window size. The width of the sliding window is adjusted according to the noise level. When the window is large, the impact of high-frequency noise can be reduced more smoothly. Traverse the logging eigenvalues, set the sliding window size to 5 data points, the window is at the center of the current point, and calculate the average value of the data points in the window as the processed eigenvalue data. For data segments with particularly obvious noise (PERM data segments), the present invention uses the median filtering method to replace the value of each data point with the median of the data points in its neighborhood, avoiding the influence of outliers on the smoothing effect. The specific formula is: Among them, the window size is k, and the median of all data points in the window is used as the value of the position; Set its window size to 5 data points, with the window centered at the current point; Next is data normalization. Since the data magnitudes of different logging characteristics vary greatly, if they are directly input into the model, the weights of some features may be too large or too small, thereby affecting the stability and convergence speed of the training. Therefore, the present invention adopts a normalization method to convert the data ranges of all features into the same scale, thereby ensuring that the feature value size does not affect the learning effect of the model. Specifically, two methods, Min-Max normalization and Z-score standardization, are adopted. For data without obvious outliers such as (AC, UPOR, SH, CAL, CN, DEN, etc.), Min-Max normalization is used to map the data to the [0,1] interval, and the calculation formula is: This method can ensure that the distribution range of the data is within a fixed interval, which is convenient for model processing; in addition, for data with extreme values, such as PERM and GR, Z-score standardization is used to further reduce the impact of extreme values ​​by adjusting the data to a distribution with a mean of 0 and a variance of 1.

6. The well logging interpretation method based on multi-task joint learning and adaptive multi-modal feature fusion according to claim 1, characterized in that: In step S3: dynamic adjacency matrix construction module: generating a spatial adjacency matrix based on the geographic spatial location of the wells; generating a geological similarity adjacency matrix based on the similarity of geological characteristics between the wells; Apply sliding depth window to generate deep sequence adjacency matrix; 301. Generate a spatial adjacency matrix. In order to capture the geographical adjacency relationship between wells, the two-dimensional spatial distance between wells is calculated based on the well location information (well ID, longitude, latitude) of the same block. The specific calculation formula is: Among them, x, y, and z are the longitude, latitude, and depth of the well location, respectively. Then, the distance threshold is set to 0.005 degrees (about 500 meters), and only the well pairs with a distance less than or equal to the threshold are retained as adjacency relationships. For the well pairs that meet the distance requirements, self-connection is avoided and the distance is inversely proportional (such as ) Calculate the edge weights and generate a spatial adjacency matrix. The weights exceeding the threshold are set to 0. At the same time, if the distance of the same well is 0, the weight will be infinite, so when the distance is set to 0, the weight is also 0; The lithology similarity is introduced into the spatial adjacency matrix. For each pair of wells in the same block, the first two digits represent the main lithology and the last two digits represent the secondary lithology according to the lithology coding. In order to speed up the calculation, the lithology similarity score can be calculated in parallel; By establishing a process pool and initializing an empty list pool_args to store all parameter pairs that need to be processed in parallel; through nested double loops, traverse all its formation lithology pairs (make sure each pair of wells is calculated only once and maintain symmetry, that is, i <= j)), extract the well ID and lithology data of each pair of wells and add them to pool_args, and then assign tasks to each process through the pool_args list to parallelly calculate the similarity scores of all formation lithologies, where the similarity scores are as follows: (1) The primary lithology and secondary lithology are the same: 1.0 (2) The primary lithology is the same, and the secondary lithology is different but similar: 0.7 (3) The primary lithology is the same, but the secondary lithology is different and dissimilar: 0.5 (4) Different primary lithology and same secondary lithology: 0.5 (5) The primary lithology is different, and the secondary lithology is different but similar: 0.3 (6) The primary lithology is different, and the secondary lithology is different and dissimilar: 0.0 Then, the lithology of multiple strata is weighted averaged according to the stratum thickness to obtain the final stratum lithology score of the well pair, and finally the calculated weighted average score is assigned to the corresponding well pair in the similarity matrix; Then, based on the adjacency matrix obtained by considering the distance, the obtained lithology similarity matrix is ​​multiplied with the previous spatial adjacency matrix to update the spatial adjacency matrix; 302. Construct a geological similarity adjacency matrix, which determines the adjacency relationship by the similarity of the 17-dimensional eigenvalues ​​of the logging data except the depth; for each feature, calculate the cosine similarity between wells: Among them, T ik and T jk are the values ​​of well i and well j on feature k respectively; weighted sum of all feature similarities is performed to obtain comprehensive similarity; the similarity threshold is set to 0.65, well pairs with similarity higher than the threshold are retained, geological similarity adjacency relationships are established, and the similarity values ​​are filled into the adjacency matrix as edge weights; 303. Generate a depth sequence adjacency matrix. To capture the dynamic association of well logging data between wells, the gradually increasing depth in the logging data is regarded as pseudo time, and a sliding depth window is used to process the depth sequence data; According to the well logging characteristic data divided and stored in S2, the depth window size is set to 10 and the sliding step size is 3; The dynamic depth warping (DTW) distance of the well pair data is calculated in each depth window in a multi-process parallel manner to measure the similarity of the depth sequence; the specific formula is: Use the DTW distance calculation function to calculate the DTW distance of all sliding window pairs between each well pair. The smaller the distance, the more similar the two windows are. Accumulate the distance of each pair of windows. If there is a window pair, return its average distance; otherwise, return infinity, and finally get the average DTW distance of all window pairs between the two wells. The threshold is dynamically adjusted according to the distribution of DTW distance, the DTW distance threshold is set to 95% of the distance, and then the adjacency matrix is ​​constructed according to the threshold; Finally, the three adjacency matrices of each block are calculated through a loop, and then the corresponding adjacency matrices of each block are spliced ​​on the diagonal, and the adjacency relationship between blocks is set to 0, and finally three global adjacency matrices are obtained.

7. The well logging interpretation method based on multi-task joint learning and adaptive multi-modal feature fusion according to claim 1, characterized in that: In step S4: Adaptive multimodal feature fusion model construction: (1) Shared feature extraction layer The well logging features and geological structure features are spliced ​​in the feature dimension to form the initial feature matrix X; the well logging features of each well have 17 dimensions, and the geological structure features have 6 dimensions, so the dimension of the spliced ​​input feature matrix X is N×23 (N is the total number of wells in all blocks); Transform X using a linear layer to map the feature dimension from 23 to a hidden layer dimension of 128; activate the transformed features using the ReLU nonlinear activation function to obtain the initial feature representation H0, whose dimension is N×128; The three adjacency matrices obtained by S3 are spatial adjacency matrices A. space , geological similarity adjacency matrix A geological and the deep sequence adjacency matrix A depth , the dimension of each adjacency matrix is ​​N×N; each adjacency matrix is ​​given a learnable weight coefficient α i (i=1,2,3), the initial value is set to Ensure that each adjacency matrix contributes equally to feature extraction in the initial stage; normalize each adjacency matrix to obtain Array The calculation formula is as follows: Among them, D i is the adjacency matrix A i The degree matrix of Use the Softmax function to normalize the weight coefficient and get the normalized weight β i : Comprehensive adjacency matrix It is obtained by weighted summing of three normalized adjacency matrices: By weighted summing of three normalized adjacency matrices, we get: Using the comprehensive adjacency matrix Perform graph convolution on the initial feature representation H0; Multiply to get the propagated features Its dimension is N×128; Then a linear transformation layer is used to map the feature dimension from 128 to 64, and a ReLU activation function is applied to generate the graph convolution layer output H1, whose dimension is N×64; In order to enhance the selectivity of the model for important features, step S4 introduces a gating mechanism; the output H1 of the graph convolutional layer passes through a linear layer and a Sigmoid activation function to generate a gating coefficient matrix G, whose dimension is also N×64 and the value range is between 0 and 1; the gated feature representation H is obtained by performing an element-wise product operation on H1 and g gated , with a dimension of N×64; this process enables the model to dynamically adjust the transfer weights of each feature channel, highlight key features, and suppress redundant information; After the gating mechanism, the model further integrates an adaptive adjacency matrix attention mechanism to optimize the weight distribution of different adjacency matrices; the model defines a new learnable parameter γ for each adjacency matrix i (i=1,2,3), the initial value is set to Through the Softmax function, γ i Normalized to weight δ i , ensure that the weight sum is 1; comprehensive adjacency matrix It is calculated by the following formula: Among them, A i is the original adjacency matrix; Then, a multi-head attention mechanism is introduced to further improve the expressiveness of features; The number of attention heads is set to 4, and each attention head processes 16 subspaces of the feature dimension; for each attention head, the model first generates a query vector, a key vector, and a value vector through linear transformation, and then calculates the attention coefficient between nodes Based on these coefficients, the features of neighboring nodes are weighted and summed to obtain the new feature representation of each node under this attention head. The outputs of all attention heads are concatenated in the feature dimension to form the final attention output H attn , whose dimension is N×64; In order to integrate the gated features H gated And the output H of the multi-head attention mechanism attn The model uses a residual connection to add the two and pass the ReLU activation function to generate the final fusion feature representation H final , whose dimension is N×64; Finally, for different prediction targets, the model designs independent multi-task output layers; for lithology prediction (classification task), a fully connected layer is used to integrate the fusion feature H final Mapped to the number of categories C, and applying the Softmax activation function, the probability distribution y of the lithology category is output lithology , whose dimension is N×C; for the permeability y prem , porosity y por and formation saturation y soi prediction (regression task), using independent fully connected layers to fuse the fusion features H final Mapped to a single scalar value, the output dimensions are all N×1.

8. The well logging interpretation method based on multi-task joint learning and adaptive multi-modal feature fusion according to claim 1, characterized in that: In step S5, multi-task joint learning: mainly through the shared feature extraction layer, the multi-task loss balance module and the multi-task prediction module to achieve the joint prediction of the three logging indicators of PEM, POR and SOI and lithology classification; the goal of this module is to effectively utilize the shared spatial features and improve the comprehensive prediction accuracy of multiple tasks, including:

501. In the shared feature extraction layer, the model sets a shared feature extraction layer for different tasks. The shared feature extraction layer is further transformed and nonlinearly activated from H final Extract higher-level feature representations, which not only include the basic features of well logging indicators such as PEM, POR, and SOI, but also incorporate the geological information required for lithology classification. By sharing the feature extraction layer, the model avoids learning feature representations for each task separately, reducing the computational burden while ensuring the consistency and complementarity of features between tasks.

502. Use adaptive adjustment of the loss weight of each task to ensure balanced optimization of each task during model training. Since the loss scale and learning rate of each task are different, directly adding their loss functions may cause the loss of a task to dominate the overall optimization process, resulting in a decrease in the prediction accuracy of the remaining tasks. To solve this problem, a weight factor α can be introduced, and the loss function L of each task task is adaptively adjusted, the formula is: L total =α1L PEM +α2L POR +α3L SOI +α4L lithology Among them, L PEM , L POR , L SOI are the mean square error losses of permeability, porosity and formation saturation prediction, L lithology is the cross entropy loss for lithology classification; weight factors α1, α2, and α3 are dynamically adjusted during model training to achieve adaptive weight allocation for each task; Adaptive weight allocation uses a dynamic weighting method to enable the model to find the best weight configuration during the continuous optimization process, ensuring the reasonable contribution of the loss of each task to the total loss and avoiding excessive impact of a certain task on the overall training. By introducing an automatic weight adjustment mechanism, α1, α2 and α3 are dynamically adjusted according to the loss changes of each task, or the weights are adjusted based on the uncertainty of the task, so that the model can achieve balanced optimization in the multi-task learning process. Finally, the model outputs to the multi-task prediction module, which contains independent output layers for multiple tasks to generate the final prediction values ​​and lithology classification results of PEM, POR and SOI. For each task, the output layer contains a fully connected layer to map the shared features to the target logging indicators. The output layer of each task is respectively calculated through the corresponding loss function L PEM , L POR , L SOI and L lithology Optimize For each task, an independent fully connected layer is designed to share the shared features H shared Map to target logging indicators or categories; the specific implementation is as follows: Permeability (PEM) prediction: A fully connected layer maps the shared features to a single scalar value and outputs the predicted value of the permeability Porosity (POR) prediction: A fully connected layer maps the shared features to a single scalar value and outputs the predicted porosity value Formation saturation (SOI) prediction: A fully connected layer maps the shared features to a single scalar value and outputs the predicted value of formation saturation Lithology classification: A fully connected layer maps the shared features to the number of categories C, and applies the Softmax activation function to output the probability distribution of lithology categories. Where CCC is the number of lithology categories; Through the multi-task prediction module, the model can independently perform regression prediction and lithology classification of PEM, POR, and SOI based on shared features, thereby achieving more accurate multi-task logging indicator prediction; the independent output layer of each task ensures that each task performs personalized prediction based on shared features, avoiding interference between tasks and improving the overall prediction performance.

9. The well logging interpretation method based on multi-task joint learning and adaptive multi-modal feature fusion according to claim 1, characterized in that: It also includes step S6 model training: Using the preprocessed and enhanced feature matrix X (dimension is N×17) obtained in steps S1 and S2, the adjacency matrix set [A space ,A geological ,A depth ] (the dimension of each adjacency matrix is ​​N×N), and the corresponding label data y lithology ,y perm ,y porosity ,y saturation , the data set is divided into training set, validation set and test set by evenly sampling the data of each block according to the block, and the ratio is set to 8:1:1; Then define the loss function of the classification task (lithology prediction) as cross entropy loss (CrossEntropyLoss), The loss function of the regression task (PEM, POR, SOI) is defined as the mean square error loss (MSELoss); The model built using S4 is trained in a sliding window loop with a size of 10 and a step size of 3, and the cumulative loss is initialized to 0; Then the feature data, three adjacency matrices, and four label data of the N wells processed each time are moved to the designated computing device to enable batch processing; Then clear the accumulated gradients in the optimizer to avoid the impact of gradient accumulation on the current step; Input the input feature X and three adjacency matrices into the model to obtain the prediction output of each task, including the prediction results of lithology classification (N, C), permeability (N, 1), and porosity (N, 1); Then the loss predicted by lithology and the predicted losses of PEM, POR, and SOI are calculated; Then the physical constraint regularization loss is calculated to ensure that the model predictions conform to actual geological laws:

601. Empirical Relationship between Porosity and Depth In reservoir geology, porosity usually decreases with increasing depth. This phenomenon is mainly due to the compaction effect of the formation. A common empirical formula used to describe the relationship between porosity and depth is as follows: φ=φ0e -c·d Where: φ represents porosity, φ0 represents the initial porosity at the shallow layer, c is the compaction coefficient (depending on the formation properties, usually between 0.0001 and 0.0005), and d represents the depth; This formula can be used for regularization constraints to make the relationship between porosity and depth predicted by the model consistent with the compaction effect; the following regularization term is added to the loss function to calculate the mean square error between the predicted porosity and the expected porosity:

602. Relationship between water saturation and depth: Water saturation S w It usually decreases with increasing depth, especially in deeper parts of the reservoir, where the water saturation tends to be lower due to capillary pressure. According to the research of JD Buckley et al. (1995), the following empirical formula can be used to describe the relationship between depth and water saturation: S w =S wc +(S wi -S wc )e -α·d Where: S w c is the residual water saturation, which represents the final water saturation value at a deeper level, S wi is the shallow water saturation value, α is the water saturation variation coefficient (usually between 0.0001 and 0.0003), and d represents the depth; This formula can be used as a regularization constraint to ensure the depth variation trend of water saturation; the regularization term is as follows, which is used to predict the mean square error between water saturation and expected water saturation:

603. Relationship between porosity and density: Formation density ρ usually increases with depth, while porosity decreases with increasing density; the commonly used relationship is as follows Where: matrix is the matrix density of the formation rock, which is usually a known value (for example, the matrix density of limestone is about 2.71 g / cm 3 );ρ fluid is the density of the fluid (such as water or oil, about 1g / cm 3 ); This formula can express the relationship between porosity and density as a regularization constraint:

604. Relationship between water saturation and relative permeability In reservoir description, water saturation S w and relative permeability k r There is a correlation between k r =k r0 ·(1-S w ) n Where: k r0 is the initial permeability when the saturation is zero, n is the permeability curve index, which usually takes a value of 2 to 4, depending on the reservoir type; This relationship can be used as a constraint in the model to ensure reasonable water saturation and permeability:

605. Complete loss function: In order to implement these regularization constraints, the above regularization terms can be added to the calculation of the loss function and the model parameters can be updated through back propagation; the overall form of the loss function is: L=L base +λ1R porosity-depth +λ2R water-saturation-depth +λ3R porosity-density +λ4R permeability-saturation L base =L lithology +L perm +L porosity +L saturation Then the gradient of the total loss L with respect to the model parameters is calculated through back propagation; Use the Adam optimizer to update the model parameters according to the calculated gradients to minimize the total loss; Add the total loss of the current batch to the total loss, and then calculate the average loss of the current epoch After each epoch, the model is evaluated using the validation set to calculate the performance indicators RMSE and accuracy of each task; Finally, according to the performance of the validation set, the loss weight α is dynamically adjusted i and the regularization weight β i .

Images