Shale reservoir cross-scale pore characterization method and system based on graph neural network

Abstract

The present application provides a shale reservoir cross-scale pore characterization method and system based on a graph neural network, belonging to the field of rock physical parameter evaluation and artificial intelligence technology, comprising: cooperatively obtaining nuclear magnetic resonance echo data and scanning electron microscope or micron CT image data of a shale sample; performing pore segmentation on the image data and constructing a pore topology graph; extracting relaxation characteristics from the echo data and performing initial inversion to obtain an initial T2 distribution; inputting the topology graph, relaxation characteristics and initial T2 distribution into a graph neural network multi-modal fusion model to generate an adaptive correction function; using the correction function to perform nonlinear reconstruction on the initial T2 distribution to obtain a corrected T2 distribution; and further outputting a full-scale pore size distribution model and multiphase fluid saturation. The present application realizes physical constraint correction of a high-resolution image on a nuclear magnetic signal through feature-level fusion, solves the multi-solution problem of ill-posed inversion of a nuclear magnetic resonance T2 spectrum, and improves the shale reservoir pore characterization accuracy.

Description

Technical Field

[0001] This invention relates to the fields of rock physical parameter evaluation and artificial intelligence technology, and in particular to a method and system for cross-scale porosity characterization of shale reservoirs based on graph neural networks. Background Technology

[0002] Shale reservoirs present significant technical challenges to accurately evaluating their physical properties due to their extreme heterogeneity, multi-scale pore systems ranging from nanometer to micrometer, and extremely complex fluid occurrence states. Currently, the industry mainly relies on three complementary detection methods to characterize the pore structure of shale reservoirs: low-field nuclear magnetic resonance (LF-NMR) technology non-destructively probes fluid information in the formation by measuring the transverse relaxation time T2 of hydrogen nuclei, indirectly inverting the pore size distribution; focused ion beam scanning electron microscopy (FIB-SEM) can provide two-dimensional or three-dimensional pore topology images with nanometer-level resolution; and micrometer-level computed tomography (Micro-CT) excels at characterizing the connectivity networks of large and medium-sized pores and fractures in three dimensions. While each of these methods has unique advantages at specific scales and dimensions, fundamental technical bottlenecks remain in cross-scale data fusion.

[0003] Nuclear magnetic resonance (NMR) T2 spectrum inversion is mathematically a typical ill-conditioned inverse problem. Because shale often contains paramagnetic minerals (such as pyrite and chlorite), these minerals generate strong internal magnetic field gradients within the pores. Combined with inherent electromagnetic noise in the instrument, this often severely obscures the weak NMR signals from tiny pores. In low signal-to-noise ratio environments, traditional mathematically regularized inversion methods (such as Tikhonov regularization and singular value decomposition) inevitably introduce multiple solutions, making it difficult to accurately recover the short T2 relaxation component signals corresponding to nanoscale pores. This leads to significant deviations in the evaluation of total porosity and oil saturation in shale reservoirs.

[0004] On the other hand, while visual imaging techniques such as SEM and Micro-CT can provide extremely precise information on the pore geometry framework, their fundamental drawback lies in the fact that these two techniques belong to "static solid phase" characterization. They can clearly observe the morphological characteristics of organic matter pores, intergranular pores, and microcracks, but lack direct information on the fluid occurrence state and dynamic flow. Correspondingly, nuclear magnetic resonance belongs to "dynamic liquid phase" characterization. It can distinguish between free oil, adsorbed oil, pore free water, and clay-bound water, but lacks clear geometric spatial constraints. How to deeply integrate high-resolution solid phase topology with low-resolution liquid phase relaxation signals has long been a challenge for the rock physics community.

[0005] Existing multi-source data fusion methods mostly remain at the stage of decision-level fusion or simple data-level splicing. For example, different logging curves are directly combined into an input vector and then fed into a multilayer perceptron or support vector machine for regression prediction, or the output results of each independent model are weighted and averaged in the final prediction stage. These loosely coupled fusion methods fail to establish the intrinsic physical relationship between different modal data at the feature level, and cannot use the structural prior of high-resolution modes to constrain and correct the inversion error of low-resolution modes. Therefore, it is difficult to fundamentally break through the information representation limit of a single mode. In view of this, there is an urgent need for an intelligent method that can achieve deep multimodal fusion and mutual correction at the feature level, so as to comprehensively and accurately reveal the true pore structure and fluid occurrence state of shale reservoirs. Summary of the Invention

[0006] To address the problems in the background technology, this invention provides a method and system for cross-scale porosity characterization of shale reservoirs based on graph neural networks. By constructing a multimodal cross-modal constraint correction modeling framework based on feature-level fusion and topological graph mapping, and utilizing the solid-state geometric features of high-resolution images to mathematically constrain and correct the liquid signal inversion of low-resolution nuclear magnetic resonance, this invention fundamentally solves the problem of multiple solutions in ill-conditioned nuclear magnetic resonance inversion, thereby comprehensively and accurately revealing the true pore structure and fluid occurrence state of shale reservoirs.

[0007] To achieve the above objectives, this invention provides a method for cross-scale porosity characterization of shale reservoirs based on graph neural networks, comprising the following steps: Collaborative acquisition of low-field nuclear magnetic resonance raw spin echo data and high-resolution scanning electron microscopy or micron-CT image data of shale samples from the same stratum; The image data is segmented to create a pore topology map structure. Relaxation features are extracted from the raw nuclear magnetic resonance spin echo data, and an initial inversion is performed to obtain the initial T2 distribution; The pore topology data, relaxation features, and initial T2 distribution are input into a multimodal fusion model of a graph neural network. Geometric features and relaxation features are aligned in the shared hidden layer space to generate an adaptive correction function. The signal distortion of the initial T2 distribution in the low signal-to-noise ratio range is nonlinearly reconstructed using the adaptive correction function to obtain the corrected T2 distribution. Based on the corrected T2 distribution, a full-scale pore size distribution model and multiphase fluid saturation are output.

[0008] As a further improvement of the present invention, when performing pore segmentation on the image data, a U-Net network with a fusion attention mechanism is adopted. The U-Net network includes an encoder path and a decoder path. The encoder path extracts multi-scale features through multi-level downsampling, and the decoder path fuses features of different scales through skip connections. The attention mechanism introduces spatial attention gating units in the skip connections. The segmentation output includes classification label maps of organic matter pores, intergranular pores, and microcracks.

[0009] As a further improvement of the present invention, the data for constructing the pore topology map structure includes: Each segmented independent pore body is abstracted as a node, and the node feature vector includes at least specific surface area, pore volume, aspect ratio, perimeter and equivalent diameter; the pore throats connecting adjacent pore bodies are abstracted as edges, and the edge feature vector includes at least pore throat radius and connectivity weight; and an adjacency matrix is ​​constructed to describe the topological connection relationship between pore bodies.

[0010] As a further improvement of the present invention, the initial inversion adopts a multi-exponential decay model based on Tikhonov regularization, and the inversion objective function is expressed as: ; The constraints are as follows: A represents the kernel matrix. Let M represent the T2 distribution vector to be determined, and let M represent the measured echo vector. Let L be the regularization parameter and L be the regularization operator.

[0011] As a further improvement of the present invention, the extraction of the relaxation features is achieved through a one-dimensional convolutional neural network, which includes three stages in sequence: convolution operation, ReLU nonlinear activation and max pooling downsampling, and outputs a compressed relaxation feature vector after multi-layer cascading.

[0012] As a further improvement of the present invention, the graph neural network multimodal fusion model adopts a graph attention network architecture, and the formula for calculating the attention coefficient is as follows: ; The node feature update rule is as follows: ; softmax represents the normalized exponential function, subscript ; j Indicates a node i All neighboring nodes j Normalization is performed so that the nodes The sum of the attention weights for all its neighbors is 1; This represents a modified linear activation function with a negative slope. ,in β The preset negative slope coefficient (0 < β<1, used to generate a non-zero output when the input is negative, in order to maintain gradient fluidity and prevent neuron deactivation; This represents a vector concatenation operation; , Represents node features, represents edge features, W, W e Represents the learnable weight matrix. Represents the attention parameter vector. This represents the transpose of the attention parameter vector, where N(i) represents the set of neighbor nodes. This represents a non-linear activation function.

[0013] As a further improvement of the present invention, feature alignment is performed in the shared hidden layer space through a contrastive learning loss function, so that the geometric features and relaxation features from the same layer and similar lithofacies are close in the hidden space, and features from different lithofacies or with significantly different pore structures are far away in the hidden space. The contrastive learning loss function adopts the InfoNCE loss, expressed as: , Among them, z geo z represents the latent space feature vector from the topological branch of the image. nmr Let represent the latent space feature vector from the NMR relaxation branch, and sim(·,·) represent the cosine similarity function. This represents the temperature hyperparameter, and V represents the number of negative samples for comparison.

[0014] As a further improvement of the present invention, generating the adaptive correction function includes: A cross-domain mapping between the T2 distribution and the pore nodes is established through a cross-attention mechanism. The embedding representation q(m) of the T2 component is used as the query vector to process the pore node embedding set output by the graph neural network. Attention-weighted aggregation is performed to obtain a comprehensive geometric feature representation. ,in, ; The adaptive correction function is expressed as follows: ; ; in, For feature dimension, This represents a correction network consisting of fully connected layers. This represents the learnable parameters of the correction network, to which a non-negative truncation operation is applied to ensure... .

[0015] As a further improvement of the present invention, the full-scale pore size distribution model adopts a splicing and fusion method based on a pore size-dependent dynamic weight function, which is expressed as follows: ; Wherein, the weight function satisfies + + =1, Obtained through differential surface relaxation rate conversion; Based on the pore type classification results from image segmentation, differentiated surface relaxation rates are assigned to organic pores, inorganic mineral pores, and microcracks. Utilizing cross attention weights β m,i Calculate the weighted equivalent surface relaxation rate ,implement Aperture transformation mapping.

[0016] This invention also provides a graph neural network-based multi-scale porosity characterization system for shale reservoirs, used to implement the aforementioned graph neural network-based multi-scale porosity characterization method for shale reservoirs, comprising: The data acquisition module is used to collaboratively acquire low-field nuclear magnetic resonance raw spin echo data and high-resolution scanning electron microscope or micron CT image data of target shale samples in the same stratum; The image segmentation module uses a U-Net network with an attention fusion mechanism to perform pore space segmentation processing on high-resolution image data, and identifies and classifies organic matter pores, intergranular pores and microcracks; The topology graph construction module constructs porous topology graph structure data containing node feature matrices and edge feature matrices based on the segmentation results; The NMR feature extraction module is used to extract features from nuclear magnetic resonance echo data to obtain relaxed feature vectors, and to perform initial T2 inversion based on Tikhonov regularization to obtain the initial T2 relaxation distribution sequence. The GNN fusion correction module is used to align geometric features and relaxation features in the shared hidden layer space through the graph neural network multimodal fusion model, generate an adaptive correction function, and perform nonlinear reconstruction and correction of signal distortion in the low signal-to-noise ratio range with the initial T2 relaxation distribution. The results output module is used to output the full-scale pore size distribution model and the quantitative evaluation results of multiphase fluid saturation.

[0017] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention, by constructing an end-to-end framework of "image segmentation-topology mapping-feature extraction-latent space alignment-adaptive correction," achieves for the first time cross-modal deep fusion of high-resolution solid-state geometric features and low-resolution liquid-phase relaxation features of shale reservoirs at the feature level. Unlike traditional data layer splicing or decision layer weighted fusion methods, this invention utilizes graph neural networks to align multimodal features in a shared hidden space, enabling high-resolution image topological information to serve as a physical prior constraint on the ill-conditioned inversion process of low-resolution NMR signals. This fundamentally solves the long-term multi-solution problem caused by noise propagation and incomplete regularization in NMR T2 spectrum inversion, significantly improving the accuracy and reliability of shale reservoir porosity characterization.

[0018] This invention employs a U-Net network with an attention-integrated mechanism to segment pores in SEM / CT images, transforming the pore network into a topological graph structure containing node feature matrices and edge feature matrices. Compared to traditional pore extraction methods based on threshold segmentation or watershed algorithms, this invention can automatically identify different pore types such as organic matter pores, intergranular pores, and microcracks. Furthermore, the graph structure naturally preserves the spatial connectivity between pores, providing a data representation that fits the physical morphology of shale pores for subsequent message passing in the graph neural network.

[0019] This invention uses a one-dimensional convolutional neural network to automatically extract relaxation features from NMR echo data, replacing the empirical parameter selection method that relies on manually designed features in traditional methods. At the same time, it uses a multi-exponential decay model based on Tikhonov regularization as the initial inversion benchmark, which takes into account both the physical interpretability of the inversion process and the end-to-end learning capability of feature extraction.

[0020] This invention employs a graph attention network to construct a multimodal fusion model. By using a contrastive learning loss function, it forces geometric and relaxed features to align in the latent space, and generates an adaptive correction function for each T2 component through a cross-attention mechanism. Compared to traditional regularized inversion methods, this mechanism can perform nonlinear reconstruction and completion of signal distortion in low signal-to-noise ratio ranges, effectively eliminating the interference of paramagnetic minerals on the T2 relaxation signal and significantly improving the characterization accuracy of nanoscale micropores.

[0021] This invention employs a stitching and fusion method based on a pore size-dependent dynamic weighting function. It adaptively fuses the nanoscale pore size distribution obtained from SEM, the micrometer-scale pore size distribution obtained from CT, and the full-scale continuous pore size distribution obtained from corrected NMR. Furthermore, it performs differentiated surface relaxation rate conversion based on the pore type classification results from image segmentation. Compared to traditional single-scale pore size characterization methods or simple weighted average fusion approaches, this invention achieves seamless continuous characterization from the nanometer to the micrometer and then to the macroscopic scale, providing reliable technical support for the accurate evaluation of multiphase fluid saturation in shale reservoirs. Attached Figure Description

[0022] Figure 1 This is a flowchart of a method for characterizing the porosity of shale reservoirs across scales based on graph neural networks, as disclosed in one embodiment of the present invention. Figure 2 This is a schematic diagram of the pore topology map construction process disclosed in one embodiment of the present invention; Figure 3 This is a network architecture diagram of a graph neural network multimodal fusion model disclosed in one embodiment of the present invention; Figure 4 This is a schematic diagram of latent space feature alignment disclosed in an embodiment of the present invention. In the figure, φ_geo(h^(L)) represents the image topology branch mapping function projecting the node embedding output by the GNN, and obtaining the sample-level feature z_geo after global mean pooling; φ_nmr(c_nmr) represents the NMR relaxation branch mapping function projecting the relaxed feature c_nmr extracted by 1D-CNN, and obtaining the sample-level feature z_nmr. Figure 5 This is a comparison diagram of the T2 distribution curve before and after correction, as disclosed in an embodiment of the present invention; Figure 6 This is a diagram illustrating the architecture of a graph neural network-based multimodal porosity characterization system for shale reservoirs, as disclosed in one embodiment of the present invention. Detailed Implementation

[0023] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0024] The present invention will now be described in further detail with reference to the accompanying drawings: like Figure 1 As shown, the shale reservoir cross-scale porosity characterization method based on graph neural networks provided by this invention includes the following steps: S1. Collaboratively acquire low-field nuclear magnetic resonance raw spin echo data and high-resolution scanning electron microscope or micron-CT image data of shale samples in the same stratum; The collaborative acquisition of multimodal raw data of the target shale geological samples includes: measuring the target shale core samples using a low-field nuclear magnetic resonance spectrometer to acquire raw spin echo attenuation signal sequences; and simultaneously performing two-dimensional / three-dimensional high-resolution imaging with focused ion beam scanning electron microscopy (FIB-SEM) and three-dimensional imaging with micron-level computed tomography (Micro-CT) on core samples from the same stratum to acquire high-resolution pore structure image data covering the nanometer to micron scale.

[0025] The acquired image data is preprocessed, including Gaussian filtering for noise reduction and normalization.

[0026] Specifically, the expression for Gaussian filtering is: I(x,y)=I 0 (x,y) G(x,y;ξ) In the formula, These are the pixel intensity values ​​of the original image. G(x,y;ξ) The standard deviation is ξ The two-dimensional Gaussian kernel function, For convolution operations, This is the filtered image.

[0027] The normalization of image pixel values ​​is represented as follows: ; In the formula, These are the original pixel values. and These are the minimum and maximum pixel values ​​in the image, respectively. The value is the normalized pixel value, and its range is [0,1].

[0028] The raw NMR spin echo data were also subjected to baseline correction and noise preprocessing to remove the effects of instrument baseline drift and high-frequency random noise.

[0029] This step establishes the correspondence between samples at the same stratum during the data acquisition stage and performs targeted preprocessing on the images and NMR signals, laying a data consistency foundation for subsequent cross-modal feature alignment and effectively avoiding spurious correlations introduced by differences in data sources.

[0030] S2. Perform aperture segmentation on the image data and construct aperture topology data; In the process of pore segmentation of image data, a U-Net network with an attention fusion mechanism is used to perform pixel-level semantic segmentation on preprocessed SEM and Micro-CT images. The U-Net network includes an encoder path and a decoder path. The encoder path achieves multi-level downsampling through alternating convolutional layers and max pooling layers, progressively extracting multi-scale features from low-level texture to high-level semantics. The decoder path achieves upsampling through transposed convolution and fuses features of different scales through skip connections. The attention mechanism introduces spatial attention gating units in the skip connections to adaptively select the encoder features most relevant to the current decoding layer, effectively suppressing the interference of background noise and irrelevant regions on segmentation accuracy. The segmentation output includes classification label maps of organic pores, intergranular pores, and microcracks.

[0031] Before extracting node and edge features, connectivity analysis is first performed on the segmented pore label map to identify independent pore bodies and their throat connections. This can be achieved using the maximum sphere method, the central axis extraction method, or the connected component labeling and throat detection method. Taking the maximum sphere method as an example, for a 3D pore image, with each voxel point in the pore space as the center of a sphere, the radius of the sphere is continuously increased until it reaches the boundary of the solid matrix, and the maximum inscribed sphere radius of each voxel point is recorded. The watershed algorithm is used to segment the pore space based on the local maximum sphere radius. Each local maximum region corresponds to an independent pore body, and the contact area between adjacent pore bodies is identified as a throat. The throat radius is taken as the average or minimum value of the maximum inscribed sphere radius of all voxel points within the contact area, and the connectivity weight is calculated by the ratio of the throat cross-sectional area to the throat length.

[0032] The construction of the pore topology graph structure data includes: abstracting each segmented independent pore body as a node, with the node feature vector containing at least specific surface area, pore volume, aspect ratio, perimeter, and equivalent diameter; abstracting the pore throats connecting adjacent pore bodies as edges, with the edge feature vector containing at least the pore throat radius and connectivity weight; and constructing an adjacency matrix to describe the topological connection relationship between pore bodies.

[0033] Furthermore, it also includes a paramagnetic mineral interference elimination step: the type and content distribution of paramagnetic minerals are incorporated as additional node features into the pore topology map, and the graph neural network learns the influence of paramagnetic minerals on the local magnetic field to adaptively cancel their shortening effect and distortion interference on the T2 relaxation signal.

[0034] The pore topology map adopts a hierarchical heterogeneous map structure: a nanoscale sub-map is constructed based on SEM images, and a micrometer-scale sub-map is constructed based on Micro-CT images. The two sub-maps interact with each other through a cross-scale aggregation operator to inject a global pooling feature summary of one scale into the node feature update process of the other scale.

[0035] In some embodiments, the segmentation network is trained using a weighted combination of the cross-entropy loss function and the Dice loss function. The cross-entropy loss function is expressed as: ; In the formula, Index for pore types (including organic matter pores, intergranular pores, microcracks, and matrix background). For pixel index, One-hot encoding for the real label, This represents the class probability predicted by the network.

[0036] The Dice loss function is expressed as: ; In the formula, Total number of categories (in this embodiment) ), For category indexing, To prevent a smoothing constant with a denominator of zero, the total loss function is: , These are the weighting coefficients for the Dice loss.

[0037] After segmentation, the output includes classification label maps containing three types of pores: organic matter pores, intergranular pores, and microcracks, providing accurate geometric boundary information for the subsequent construction of pore topology maps.

[0038] like Figure 2 As shown, based on the segmentation results, the geometric structure of the shale pore network is transformed into a mathematical graph. In this graph structure, each independent pore volume is abstracted as a node, and the throats connecting adjacent pore volumes are abstracted as edges.

[0039] For each node, a multidimensional feature vector representing its geometric shape is extracted. The feature vector of each node is defined as follows: ; In the formula, For the first Specific surface area of ​​each porous body For pore volume, The aspect ratio (the ratio of the major axis to the minor axis). The perimeter of the pore boundary. The equivalent circle diameter is used; the resulting node feature matrix is ​​used to provide the graph neural network with solid-state geometric property information for each pore body.

[0040] For each edge, a feature vector representing its connectivity is extracted. The feature vector of each edge is defined as follows: ; In the formula, To connect porous bodies and The radius of the throat of the orifice, Connectivity weights are defined as follows: ,in To connect porous bodies and The length of the pore throat is defined as being physically proportional to the fluid conduction capacity of the pore throat channel; the resulting edge feature matrix is ​​used to provide information on the spatial connectivity and fluid transport capacity between pores to the graph neural network.

[0041] The node feature matrix X has a dimension of N×F, and the edge feature matrix E has a dimension of M×G. An adjacency matrix A is constructed to describe the topological connections between porous structures. Its dimension is N×N. That is, the adjacency matrix of the graph is defined as follows: (If porous body) and (There is a throat connection between them). (other); The node feature matrix and edge feature matrix constructed in this way are respectively expressed as: ; ; In the formula, N is the total number of nodes (pore bodies) in the graph, F is the feature dimension of each node (F=5 in this embodiment), M is the total number of edges (throat) in the graph, and G is the feature dimension of each edge (G=2 in this embodiment).

[0042] When constructing the pore topology map, the nanoscale pore network obtained by SEM and the microscale pore network obtained by Micro-CT are mapped independently, forming nanoscale and microscale pore topology sub-maps. Considering the order-of-magnitude difference in the observation scale of SEM images (field of view typically a few micrometers to tens of micrometers) and Micro-CT images (field of view typically hundreds of micrometers to millimeters) in physical space, rigid physical space edge connectivity is not mandatory between the two scale sub-maps. Instead, a cross-scale aggregation operator is used to perform statistically significant high-dimensional feature interaction: the features of all nodes in the microscale sub-map are global pooled (such as mean pooling or attention pooling) to generate a microscale global feature summary vector. This global feature summary is injected as additional contextual information into the GNN feature update process of each node in the nanoscale sub-map, and vice versa. This overcomes the engineering bottleneck of achieving absolute alignment of microscopic coordinates during sample preparation and scanning. The resulting hierarchical heterogeneous graph structure maintains the topological independence of each scale sub-map while achieving effective fusion of cross-scale statistical features.

[0043] This step transforms the shale pore network into a topological graph structure containing rich geometric information, and uses a hierarchical heterogeneous graph to achieve cross-scale feature interaction between nanoscale and microscale subgraphs. This overcomes the engineering bottleneck of the difficulty in absolutely aligning the coordinates of physical samples. At the same time, paramagnetic minerals are incorporated as additional node features into the graph structure, providing structured prior information for subsequent graph neural networks to eliminate mineral interference.

[0044] S3. Extract relaxation features from the raw spin echo data of nuclear magnetic resonance and perform initial inversion to obtain the initial T2 distribution; The extraction of relaxation features is achieved through a one-dimensional convolutional neural network (1D-CNN). This network consists of three stages: convolution operation, ReLU nonlinear activation, and max pooling downsampling. After multi-layer cascading, it outputs a compressed relaxation feature vector. The one-dimensional convolutional neural network adopts a cascaded structure of 3 to 5 convolutional layers, with the specific number of layers determined based on the length and signal complexity of the NMR echo data. For conventional acquisition parameters with the number of echo points N between 1024 and 4096, a 4-layer convolutional structure is preferred; for fast acquisition modes with fewer echo points (N < 512), a 3-layer convolutional structure is preferred; and for high-resolution acquisition modes with more echo points (N > 8192), a 5-layer convolutional structure is preferred.

[0045] The initial inversion uses a Tikhonov-regularized multi-exponential decay model, and the inversion objective function is expressed as: ; In the formula, This represents the vector that minimizes the objective function. ; Denotes the square of the L2 norm; the constraint is: A represents the kernel matrix, whose elements are... , For the first Each echo time point; Let T2 be the distribution vector to be determined. ,in, M represents the total number of discrete components in the T2 spectrum, where T denotes vector transpose and M represents the measured echo vector. , where N is the total number of callback data points; This is the regularization parameter, used to balance the weights of the data fitting term and the regularization term; For regularization operators, the identity matrix or first-order difference matrix is ​​typically chosen; the solution process applies... The non-negativity constraint ensures the physical rationality of the amplitude of each component of the T2 distribution, and the resulting initial T2 relaxation distribution sequence is used as one of the inputs for subsequent GNN fusion correction.

[0046] In some embodiments, during the convolution operation stage, multiple one-dimensional convolution kernels are used to perform sliding window inner product operations along the time dimension of the echo signal, converting the echo attenuation mode within each local time window into the corresponding feature response value. The convolution operation of 1D-CNN is represented as follows: ; In the formula, x ( ) is the first of the input one-dimensional NMR echo signal sequences. One sampling point, =0,1,…, -1; w (u) The first convolution kernel u One weight parameter, u =0,1,…, U -1; U The kernel size is [size]. For the bias term, y( The output is the convolutional kernel; in this embodiment, the kernel size of each convolutional layer is... U The (sliding window length) value ranges from 3 to 11, preferably 5 or 7, for smaller convolutional kernels ( U =3) is suitable for extracting local detail features, and larger convolutional kernels ( U =7 or 11) is suitable for capturing the medium- to long-term attenuation trend of echo signals. Each layer can use the same size convolution kernel, or a layer-by-layer decreasing strategy can be used (e.g., the first layer...). U =7, intermediate layer U =5, top layer U =3).

[0047] In the nonlinear activation phase, the ReLU activation function is used. A nonlinear transformation is applied to the convolution output to enhance the network's ability to express nonlinear features in the relaxed decay signal; In the pooling downsampling stage, the convolutional output undergoes a non-linear transformation using the ReLU activation function, followed by downsampling (dimensionality reduction and compression) of the activated feature map through a max pooling layer. This reduces the feature dimension and extracts key relaxation decay features while preserving key relaxation decay pattern characteristics. After multiple convolution-activation-pooling operations, the compressed relaxation feature vector is finally obtained. Where D is the feature vector dimension; specifically, each convolutional layer is followed by a max pooling layer with a pooling window size of 2 and a stride of 2. That is, the feature map length is halved after each pooling layer. After 3-5 convolutional-pooling cascades, the echo signal length is compressed from the original N points to Point, where L is the number of pooling layers.

[0048] Simultaneously, an initial inversion of the raw NMR echo data was performed using a multi-exponential decay model based on Tikhonov regularization. The NMR multi-exponential decay model is expressed as: ; In the formula, M(t) is the amplitude of the echo signal measured at time t. For the first The signal amplitude corresponding to each T2 component (i.e., the i-th component of the T2 distribution). For the preset first Each relaxation time component value This represents the total number of discrete components in the T2 spectrum (usually set to 64 to 256 log-spaced components). This includes error terms that incorporate instrument noise and environmental interference.

[0049] This step automatically learns relaxation features from the original echo signal using 1D-CNN, avoiding the limitations of traditional methods that rely on manually designed features. At the same time, it retains the physical inversion based on Tikhonov regularization as the initial T2 distribution benchmark, taking into account both the flexibility of end-to-end learning and the prior constraints of the rock physics model, providing input with both physical rationality and feature representation capability for subsequent multimodal fusion correction.

[0050] S4. Input the pore topology data, relaxation features, and initial T2 distribution into the multimodal fusion model of the neural network, align the geometric features and relaxation features in the shared hidden layer space, and generate an adaptive correction function. The graph neural network multimodal fusion model employs a graph attention network architecture, and the attention coefficient is calculated using the following formula: ; In the formula, the node feature update rule is: Represents the normalized exponential function, subscript Indicates a node All neighboring nodes Normalization is performed so that the nodes The sum of the attention weights for all its neighbors is 1; This represents a modified linear activation function with a negative slope. ,in The preset negative slope coefficient (0 < <1, used to generate a non-zero output when the input is negative, in order to maintain gradient fluidity and prevent neuron deactivation; This represents a vector concatenation operation; , Representing node characteristics, Representing edge features, , Represents the learnable weight matrix. Represents the attention parameter vector. Represents the set of neighboring nodes. This represents a non-linear activation function.

[0051] In the shared hidden layer space, feature alignment is performed through a contrastive learning loss function, which makes the geometric and relaxation features from the same layer and similar lithofacies close in the hidden space, while features from different lithofacies or with significantly different pore structures are far away in the hidden space. To force the latent space representations of the two branches to be semantically aligned, the contrastive learning loss function adopts the InfoNCE loss, denoted as: ; In the formula, z geo z represents the latent space feature vector from the topological branch of the image. nmr Let sim(·,·) represent the latent space feature vector from the NMR relaxation branch, τ represent the temperature hyperparameter, V represent the number of contrast negative samples, and z represent the number of contrast negative samples. v Let V represent the latent space vector (v=1,2,3…,V) of the k-th randomly sampled negative sample. The summation of the denominator iterates through the positive samples and all V negative samples to form a set of V+1 samples. The principle for constructing positive sample pairs is: paired z geo and z nmr Samples from the same stratigraphic level, similar lithofacies, and possessing statistical isomorphism represent a statistically significant correspondence between the macroscopic average relaxation characteristics within the NMR measurement volume region and the statistical geometric characteristics of the image sampling slices within that region. Negative sample pairs consist of feature sets from stratigraphic intervals with significantly different lithofacies or pore structures. This sample pairing strategy, based on statistical representativeness rather than a one-to-one correspondence with spatial location, aligns with the physical scale differences between NMR volumetric measurements (centimeter-scale) and SEM / CT local imaging (micrometer- to millimeter-scale).

[0052] In some embodiments, a two-stage training strategy is adopted: first, the U-Net image segmentation network is trained independently until convergence, then the segmentation network parameters are frozen, and the pore topology map is constructed using its output segmentation results, and then the GNN multimodal feature fusion model in this step is trained.

[0053] like Figure 3 As shown, a multimodal feature fusion model is constructed using a graph attention network (GAT) to deeply fuse the obtained pore topology graph structure data with the NMR relaxation feature vector and the initial T2 distribution in a shared hidden layer space.

[0054] Graph attention networks aggregate features on the porous topology graph through a message-passing mechanism. In the l-th layer message passing, the attention coefficient is calculated using the following formula: ; In the formula, and The first Layer nodes and nodes eigenvectors, For the learnable linear transformation matrix of node features, The original feature dimension of the node. The transformed node feature dimension (hidden layer dimension); Let be the learnable linear transformation matrix for edge features. The original feature dimension of the edge. The dimension of the transformed edge features; For connecting nodes and The edge feature vectors, Let || be the attention parameter vector, and || be the feature concatenation operation. For nodes In the set of neighboring nodes in the graph, LeakyReLU is a modified linear activation function with a negative slope; exp is the natural exponential function. Represents the real number field.

[0055] Based on the attention coefficient, the node feature update rule is as follows: ; In the formula, ∑ is an ELU or ReLU nonlinear activation function used to perform a nonlinear transformation on the aggregated features; ∑ is the summation symbol. Indicates a node All neighboring nodes Perform a traversal and summation; is a learnable linear transformation matrix for node features; For nodes In the The updated feature vectors are processed through L layers of message passing; each node then obtains a high-level embedding representation that incorporates its multi-hop neighborhood topology information. .

[0056] like Figure 4 As shown, to achieve the alignment of multimodal features in the latent space, latent space mapping functions are defined for the image topology branch and the NMR relaxation branch, respectively. The mapping function for the image topology branch is: ; The mapping function for the NMR relaxation branch is: ; In the formula, W geo b geo Wnmr b nmr These are the learnable parameters of the two mapping functions. This represents the relaxed feature vector extracted from raw NMR echo data by a one-dimensional convolutional neural network. The latent space mapping function representing the topological branch of an image. The latent space mapping function represents the NMR relaxation branch; the node-level latent space vector is obtained by projecting it onto each pore node. The projection results of all nodes are then subjected to global mean pooling. Aggregate into sample-level feature vectors. The pooled sample-level representation; z nmr is the sample-level latent space feature vector of the NMR relaxation branch.

[0057] In this invention, the total loss function for training the image segmentation network and the GNN multimodal feature fusion model is: ; In the formula, L recon L is the T2 distribution reconstruction loss (measuring the error between the corrected T2 distribution and the reference standard). NCE To compare the learning alignment loss, L reg For regularization terms of model parameters, α and β These are the weighting coefficients for the corresponding loss terms.

[0058] L recon The reference standard T2 distribution is obtained through the following methods: For the training sample set, high-quality NMR measurement results obtained under high signal-to-noise ratio conditions (such as long-term superposition acquisition or high magnetic field intensity environment) are selected as the reference benchmark; In the case of lack of high-quality measured data, the reference T2 distribution can also be obtained by back-calculating the pore size distribution measured by high-pressure mercury intrusion experiments on cores of the same stratum, combined with the calibrated surface relaxation parameters.

[0059] This invention utilizes the multimodal association knowledge learned by a graph attention network in a shared latent space to generate an adaptive correction function for correcting the initial T2 inversion results. It is particularly noteworthy that the node embedding set output by the GNN { } subscript For pore node index ( =1,...,N), while the index m of the T2 distribution is the index of the relaxation time component (m=1,...,N). The two domains have completely different dimensions and physical meanings. Therefore, a cross-attention mechanism is introduced to establish a cross-domain mapping relationship between them.

[0060] First, the m-th relaxation time component in the initial T2 distribution is used to generate a query vector q(m) through an embedding layer. Then, q(m) is embedded with the pore nodes through a linear projection layer. Mapped to all dimensions After finding the same feature space, using q(m) as the query, embedding all pore nodes. Using the key value, a comprehensive geometric feature representation matching the T2 component is generated through scaled dot product attention-weighted aggregation: ; ; In the formula, q(m) is the query embedding vector of the m-th T2 component. For the final layer of GNN Embedding vectors of each pore node Let be the cross-attention weight of the m-th T2 component to the i-th pore node. is the feature dimension (used for scaling to prevent excessively large inner product values), and N is the total number of pore nodes. Through this mechanism, each T2 component can adaptively search for and aggregate the geometric features most relevant to its relaxation time from the global pore topology information.

[0061] Based on the comprehensive geometric feature representation, the correction function is expressed as: ; ; Among them, f initial (m) represents the amplitude of the m-th T2 component obtained from the initial inversion, f corrected (m) represents the amplitude of the m-th T2 component after correction, and Δ(m) represents the adaptive correction amount of the m-th component. The correction network consists of fully connected layers. The comprehensive geometric feature representation obtained by cross-attention aggregation, z nmr Let θ be the latent space eigenvector of the NMR branch. c To correct the learnable parameters of the network, m=1,2,..., The correction results must meet the following requirements. The non-negativity constraint is implemented by applying a ReLU truncation operation at the output of the correction network to ensure the physical rationality of the amplitude of each component of the T2 distribution.

[0062] This step aligns solid-state geometric features and liquid-phase relaxation features in a shared hidden space using a graph attention network. It utilizes contrastive learning to force positive sample pairs to be close and negative sample pairs to be far apart, achieving cross-modal physical constraints at the feature level for the first time. Furthermore, it establishes a cross-domain mapping between the T2 component and the pore nodes through a cross-attention mechanism, generating a component-wise adaptive correction function, which provides a mathematical basis for the accurate reconstruction of low signal-to-noise ratio signals.

[0063] S5. Use an adaptive correction function to nonlinearly reconstruct the signal distortion of the initial T2 distribution in the low signal-to-noise ratio range to obtain the corrected T2 distribution. Among them, such as Figure 5 As shown, in the low signal-to-noise ratio range (corresponding to the short T2 relaxation components of nanoscale micropores), the initial mathematical inversion often results in severe signal loss and amplitude distortion due to noise overload and the ambiguity of the ill-conditioned inverse problem. The adaptive correction function of this invention, through a cross-attention mechanism, enables each T2 component to search and aggregate high-resolution geometric prior knowledge most relevant to its relaxation time from the global pore topology information provided by SEM / CT. This means that these nanoscale pores do exist in the solid-state image and have definite geometric parameters, adaptively completing lost short T2 signal segments or correcting distorted T2 distribution curve segments during the inversion process. This correction process is a nonlinear reconstruction performed component-by-component of T2, with the correction amount for each component determined jointly by the global topological features aggregated through cross-attention and the NMR latent space information.

[0064] Furthermore, a judgment method based on a signal-to-noise ratio (SNR) threshold or a judgment method based on a T2 time threshold can be used to quantize and determine the low SNR interval. Taking the judgment method based on the SNR threshold as an example, the global SNR of the raw NMR echo data is calculated and defined as the ratio of the echo first-point amplitude to the noise standard deviation. The SNR threshold is then set. th The value ranges from 10 to 50, with a preferred value of 20. When the global SNR is lower than this threshold, full-spectrum correction is enabled; when the global SNR is not lower than this threshold, further component-by-component evaluation is performed: for each T2 component, when its amplitude is lower than 3-5 times the noise level, the component is determined to be in the low signal-to-noise ratio range.

[0065] This step uses an adaptive correction function to nonlinearly reconstruct and complete the short T2 components in the initial T2 distribution that are affected by noise submersion and ill-conditioned inversion ambiguity. This allows the clearly existing nanoscale pore geometry information in high-resolution SEM / CT images to constrain the correction direction of the NMR signal, fundamentally solving the problem of severe signal loss or distortion in the low signal-to-noise ratio range of traditional methods.

[0066] S6. Output a full-scale pore size distribution model and multiphase fluid saturation based on the corrected T2 distribution.

[0067] Based on the obtained corrected T2 relaxation distribution and combined with the pore geometry information directly provided by SEM and CT images, a full-scale pore size distribution model is generated using a stitching and fusion method based on a pore size-dependent dynamic weighting function. ; In the formula, w nano (r), w micro (r) and w nmr (r) represents the adaptive fusion weighting function that dynamically changes with the pore radius r, assigning PSD in the nanometer range. SEM (r) Higher weighting to leverage the direct geometric resolution of SEM, assigning PSD in the micrometer range. CT (r) Higher weighting to leverage the 3D volumetric statistical advantages of CT, utilizing PSD in overlapping regions. NMR_corrected (r) serves as a fluid volume conservation benchmark for smooth transition, and for any aperture r, the weighting function satisfies + + =1, This serves as a benchmark for the continuous pore size distribution across the entire scale from nanoscale to macroscale, obtained through differential surface relaxation rate conversion based on the corrected T2 distribution (where the short T2 interval at the nanoscale has been signal-completed and reconstructed by the GNN cross-attention correction function). Based on the pore type classification results from image segmentation, differentiated surface relaxation rates are assigned to organic pores, inorganic mineral pores, and microcracks. Utilizing cross attention weights β m,i Calculate the weighted equivalent surface relaxation rate ,implement Aperture transformation mapping.

[0068] In some embodiments, the weighting function adopts a smooth transition form based on the sigmoid function: , , ,in The sigmoid function is used, where r1 and r2 are the transition center apertures between the nanometer and micrometer scales, respectively (r1 = 100 nm, r2 = 1 μm in this embodiment). δ1 and δ2 are scale parameters controlling the transition band width. Logarithmic scaling ensures a uniform and smooth transition across apertures spanning several orders of magnitude. The values ​​of scale parameters δ1 and δ2 should ensure that the transition band width is much smaller than the logarithmic distance between r1 and r2, such that w nano (r) and w micro The effective transition intervals of (r) do not significantly overlap. In actual calculations, for w... nmr (r) Apply nonnegative truncation wnmr (r)=max(0,1-w nano (r)-w micro (r)) and renormalize the three weight functions to ensure that the sum of the weights at any aperture is strictly equal to 1. It should be noted that the PSD of each component is expressed as a normalized probability density function, meaning its integral over the entire aperture domain is 1. Since each component's PSD is independently normalized and the dynamic weights are functions of the aperture, the weighted fused PSD... total The integral value of (r) is not strictly 1, therefore renormalization is required: Then multiply by the total porosity φ_NMR obtained from the total amplitude calibration of the original NMR echo signal (independent of the T2 inversion and correction process), i.e. To restore the absolute pore volume distribution and ensure It satisfies the physical volume conservation constraint.

[0069] Differential surface relaxation rate conversion fully utilizes the pore type classification results from image segmentation. In shale, the organic matter pore surface is composed of oleophilic kerogen, and its surface relaxation rate... Typically small; the surface of inorganic mineral pores consists of hydrophilic clay minerals or quartz particles, with a surface relaxation rate of... The values ​​are typically large, sometimes differing by an order of magnitude. The surface relaxation rate parameter set can be obtained based on NMR-mercury intrusion porosimetry experiments using core samples from the same block, or it can be jointly optimized during model training as learnable parameters. Therefore, for each component in the corrected T2 distribution, the attention weights already calculated using the cross-attention mechanism are... β m,i The expected proportion of different pore types in the T2 component is obtained by summing the weights of nodes belonging to the same pore type k. (i.e., the proportion of different pore types in each T2 component), and the weighted equivalent surface relaxation rate is calculated accordingly. (k∈{org,inorg,frac}), execute Differential aperture transformation mapping (where each) The parameters have incorporated the pore geometry factor, thus obtaining a full-scale continuous pore size distribution that better reflects the actual shale rock physics.

[0070] Based on the full-scale pore size distribution model and the corrected T2 distribution, the multiphase fluid saturation was further calculated. Oil saturation and bound water saturation are expressed as follows: ; ; in, For the magnitude of the m-th component of the corrected T² distribution, the summation symbol in the above formula is... and These represent the conditions T2(m) > T for all cases where T2(m) > T. 2c_ow and T2(m) <T 2c_bw Summing the T2 component index m, T 2c_ow The T2 cutoff value (the relaxation time threshold that distinguishes between oil-bearing and water-bearing pores) is the oil-water cutoff value. 2c_bw The T2 cutoff value for bound water (the relaxation time threshold that distinguishes between bound water and movable water) must be set to meet the following conditions. The constraints are set to ensure the rationality of the physical meaning. The cutoff value is set based on the oil-bearing signal range determined by the T2 distribution combined with two-dimensional nuclear magnetic resonance experiments or core calibration. o S represents the oil saturation. wirr For bound water saturation. Free water saturation is determined by... Obtained through calculation.

[0071] Finally, the system outputs a full-scale high-precision pore size distribution model after feature-level cross-modal fusion correction, as well as quantitative evaluation results of multiphase fluid saturation including oil saturation, bound water saturation, and free water saturation, providing a scientific basis for the fine evaluation of shale reservoirs and the design of development plans.

[0072] This step employs a pore size-dependent dynamic weighting fusion strategy, assigning higher weights to SEM pore size distribution in the nanoscale range, higher weights to CT pore size distribution in the micrometer range, and using the corrected NMR pore size distribution as a volume conservation benchmark for a smooth transition across the entire scale range. This achieves seamless and continuous characterization from the nanoscale to the micrometer and then to the macroscale. Combined with differentiated surface relaxation rate conversion, the pore size conversion is made more consistent with the physical reality of different pore types in shale, providing a high-precision quantitative basis for evaluating the oil-bearing properties of shale reservoirs and designing development plans.

[0073] like Figure 6 As shown, the present invention also provides a shale reservoir cross-scale porosity characterization system based on graph neural networks, comprising: The data acquisition module is used to collaboratively acquire low-field nuclear magnetic resonance raw spin echo data and high-resolution scanning electron microscope or micron-CT image data of target shale samples in the same stratum; specifically, the data acquisition module includes an NMR data acquisition unit and an image data acquisition unit, which are respectively connected to the data output ports of the low-field nuclear magnetic resonance spectrometer and the scanning electron microscope / CT equipment.

[0074] The image segmentation module utilizes a U-Net network with an attention fusion mechanism to segment high-resolution image data into pore spaces, identifying and classifying organic pores, intergranular pores, and microcracks; this module is equipped with a pre-trained image segmentation neural network model.

[0075] The topology graph construction module constructs a pore topology graph structure data containing node feature matrices and edge feature matrices based on the segmentation results. This module implements the mapping transformation from pore bodies to nodes and from pore throats to edges, as well as the calculation of the adjacency matrix.

[0076] The NMR feature extraction module is used to extract features from nuclear magnetic resonance echo data through a one-dimensional convolutional neural network to obtain relaxed feature vectors, and to perform initial T2 inversion based on Tikhonov regularization to obtain the initial T2 relaxation distribution sequence. The GNN fusion correction module is used to align solid geometric features and relaxation features in the shared hidden layer space through a graph neural network multimodal fusion model, generate an adaptive correction function, and perform nonlinear reconstruction and correction of signal distortion in the low signal-to-noise ratio range with the initial T2 relaxation distribution. This module is the core processing unit of this system and realizes feature-level deep fusion of multimodal data.

[0077] The results output module is used to output a full-scale pore size distribution model and quantitative evaluation results of multiphase fluid saturation, including oil saturation, bound water saturation, and free water saturation, based on the corrected T2 relaxation distribution and using a splicing and fusion method based on a pore size-dependent dynamic weighting function. This module supports data visualization and report export functions.

[0078] Example: In a certain shale oil and gas exploration block, pore structure analysis and oil-bearing potential evaluation of a key exploration well's shale core are required. The shale reservoir in this block has the following characteristics: well-developed organic matter pores (nanoscale), high clay mineral content (including paramagnetic pyrite), and strong heterogeneity in pore structure. The method of this invention is used to characterize the pores of this core across multiple scales and modes, specifically including: Step 1.1 Data Acquisition: A 1-inch diameter, 3-cm length shale core plunger sample was drilled from a depth of 3300m in the target well. The core was divided into two parts: one part was used for low-field nuclear magnetic resonance (NMR) measurements, and the other part was used for SEM and Micro-CT imaging. A 2MHz low-field NMR spectrometer was used, with echo interval TE=0.2ms, echo number NE=2048, waiting time TW=6s, and 128 superimposed scans to acquire the raw spin echo attenuation signal sequence. Simultaneously, blank noise signals without samples were acquired for subsequent signal-to-noise ratio (SNR) evaluation.

[0079] Step 1.2 Image Acquisition: A small plunger with a diameter of 5 mm and a height of 5 mm was cut from the core sample of the same stratum and subjected to Micro-CT scanning (resolution 1 μm / voxel, scanning field of view 2 mm × 2 mm × 2 mm) and FIB-SEM scanning (resolution 10 nm / pixel, scanning field of view 10 μm × 10 μm × 10 μm).

[0080] Step 1.3 Preprocessing: Gaussian filtering is performed on the Micro-CT and SEM images for noise reduction. ξ =1.5) and grayscale normalization; baseline correction was performed on the NMR echo data to subtract instrument zero drift.

[0081] Step 2.1 Image segmentation: A pre-trained attention U-Net network was used to segment the SEM and Micro-CT images pixel by pixel. The segmentation results identified 342 organic pores, 156 intergranular pores, and 23 microcracks, and output classification label maps for the three types of pores.

[0082] Step 2.2 Connectivity Analysis: The maximum sphere method was used to extract pore connectivity relationships from the 3D Micro-CT images. Twelve connected pore clusters were identified, among which the largest cluster contained 187 pore bodies; Pore throat identification results: a total of 426 pore throats were identified, with pore throat radii ranging from 0.5 to 15 μm.

[0083] Step 2.3 Topology Graph Construction: Using each independent pore volume as a node (521 nodes in total), calculate the specific surface area, pore volume, aspect ratio, perimeter, and equivalent diameter of each node, forming a 521×5 node feature matrix X; using each pore throat as an edge (426 edges in total), calculate the throat radius and connectivity weight, forming a 426×2 edge feature matrix E; construct a 521×521 adjacency matrix A, where A... ij =1 indicates that pores i and j are connected.

[0084] Step 2.4 Cross-scale subgraphs: A nanoscale subgraph (234 nodes, 189 edges) is constructed based on SEM images, and a microscale subgraph (287 nodes, 237 edges) is constructed based on Micro-CT. Feature interaction between the two subgraphs is achieved through a cross-scale aggregation operator.

[0085] Step 3.1 Feature Extraction: Input the 2048 echo point data into a 4-layer 1D-CNN (kernel configuration: [24,48,96,192], kernel size: [7,5,5,5]). After convolution, ReLU activation and max pooling, output a 128-dimensional relaxed feature vector c. nmr .

[0086] Step 3.2 Initial T2 Inversion: Using the Tikhonov regularization method, set the number of T2 components. =64 (logarithmically evenly spaced distribution, ranging from 0.01ms to 1000ms), and the regularization parameter λ was determined to be 0.03 using the L-curve criterion. The initial T2 distribution with 64 components was obtained through inversion. Observations revealed that the amplitudes of short T2 components (corresponding to nanoscale organic pores) with T2 < 1ms were significantly low, and some components showed negative or zero values, reflecting a serious problem of multiple solutions in the inversion.

[0087] Step 4.1 GNN Fusion Model: Combine the topology graph (node ​​features X, edge features E, adjacency matrix A) and relaxation features c nmr A 4-layer graph attention network (GAT) is input with an initial T2 distribution (64 dimensions) and a 128-dimensional (128-dimensional) T2 distribution. GAT hyperparameters: hidden layer dimensions [64, 128, 256, 128], number of multi-head attention heads 4, Dropout=0.2.

[0088] Step 4.2 Latent Space Alignment: This is achieved through contrastive learning loss (InfoNCE, temperature). =0.07, negative sample number K=32) latent space features z of forced image branching geo Latent space features z of NMR branches nmr Alignment. After 200 rounds of training, the loss was verified to be stable at 0.23, indicating that the model converged.

[0089] Step 4.3 Adaptive Correction Function Generation: Using a cross-attention mechanism, the attention weights are calculated with the embedding representation q(m) of the T2 component as the query and the embeddings of the 521 pore nodes output by GAT as the keys. β m,i The aggregation yields a comprehensive geometric feature representation for each T2 component. .based on and z nmr The correction network (3 fully connected layers, 256 hidden layers) outputs the correction amount Δ(m) for each component.

[0090] Step 5.1 Low SNR Range Determination: The global SNR calculation value is 18 (threshold SNR) th =20), the entire spectrum is determined to be in a low signal-to-noise ratio state, and the correction mode is activated. Further division based on the T2 time threshold: 15 components with T2 < 0.5 ms are judged as the strong low signal-to-noise ratio interval, 32 components with T2 between 0.5 and 10 ms are judged as the medium low signal-to-noise ratio interval, and 17 components with T2 > 10 ms are judged as the high signal-to-noise ratio interval.

[0091] Step 5.2 Component-by-component correction: Apply λ to the 15 components where T2 < 0.5 ms corr A strong correction of 1.0 is applied to the component of T2 in the 0.5–10 ms range with λ. corr A gradual weakening correction of 0.3-0.8 is applied to the components with T2>10ms, and λ is applied to them. corr =0 (uncorrected). After correction, the previously missing or distorted short T2 components were successfully reconstructed, and the T2 distribution curve recovered its bimodal characteristics: the first peak was located at 0.3 ms (organic pores), and the second peak was located at 15 ms (intergranular pores and microcracks).

[0092] Step 6.1 Fusion of full-scale pore size distribution: A pore size-dependent dynamic weighting function was used for fusion. SEM pore size distribution (1-100 nm) had a weight of 0.7-0.9 in the nanometer range, CT pore size distribution (1-50 μm) had a weight of 0.6-0.8 in the micrometer range, and the corrected NMR pore size distribution (0.01-1000 nm+ μm continuous) served as a volume conservation baseline (weight 0.2-0.5) for smooth transition in the overlapping range. After fusion, PSD_total(r) exhibited a three-peak characteristic: a main peak at 30 nm (organic matter pores), a secondary peak at 2 μm (intergranular pores), and a small peak at 20 μm (microcracks). Multiplying by the total porosity... After reaching 6.8%, the absolute pore volume distribution was obtained.

[0093] Step 6.2 Differentiated Surface Relaxation Rate Transformation: Based on the image segmentation results, assign... =1.5μm / s (organic pores) =5.0μm / s (inorganic mineral pores) =8.0 μm / s (microcrack). Utilizing cross-attention weights β m,i Calculate the weighted equivalent surface relaxation rate of each T2 component and perform pore size transformation mapping to make the pore size distribution more consistent with the physical reality of shale.

[0094] Step 6.3 Multiphase fluid saturation calculation: Based on the corrected T2 distribution, set the oil-water cutoff value T. 2c_ow =12ms (calibrated by two-dimensional NMR experiment), bound water cutoff value T 2c_bw =2ms. The calculated oil saturation S o =43.2%, bound water saturation S wirr =31.5%, free water saturation S wf =25.3%.

[0095] Advantages of this invention: This invention achieves feature-level deep fusion of SEM / CT solid-state geometric features and NMR liquid-phase relaxation features in a shared hidden space using graph neural networks, overcoming the limitations of traditional decision-level fusion or data-level stitching methods. Traditional regularized inversion methods (such as Tikhonov regularization and singular value decomposition) inevitably introduce ambiguity in low signal-to-noise ratio environments, particularly leading to severe loss of short T2 relaxation component signals corresponding to nanoscale pores. This invention utilizes high-resolution image topological information as a physical prior to constrain the ill-conditioned inversion process of low-resolution NMR signals, enabling high-resolution solid-state geometric features to adaptively correct low-resolution liquid-phase relaxation features. This fundamentally solves the long-standing ambiguity problem in NMR T2 spectrum inversion, significantly improving the accuracy and reliability of shale reservoir porosity characterization.

[0096] This invention transforms the complex, irregular topology of shale pore networks into a mathematical graph structure. It utilizes the message-passing mechanism of Graph Attention Networks (GAT) to naturally process non-Euclidean spatial data, overcoming the inherent limitation of traditional convolutional neural networks, which can only handle regular grid data. By abstracting each independent pore volume as a node (containing geometric parameters such as specific surface area, pore volume, aspect ratio, perimeter, and equivalent diameter) and the pore throat as an edge (containing pore throat radius and connectivity weight), and constructing an adjacency matrix to describe the topological connections between pores, this invention achieves a precise characterization of the connectivity and fluid transport capacity of shale pore networks, better reflecting the actual physical morphology of shale reservoir pore structures.

[0097] This invention introduces a cross-attention mechanism to establish a cross-domain mapping relationship between the T2 distribution and pore nodes, generating a component-wise adaptive correction function. This function can perform nonlinear reconstruction and completion of signal distortion in the low signal-to-noise ratio range (corresponding to the short T2 relaxation component of nanoscale micropores). Simultaneously, by incorporating the type and content distribution of paramagnetic minerals as additional node features into the pore topology graph, the graph neural network learns the influence of paramagnetic minerals on the local magnetic field during message passing, adaptively offsetting the shortening effect and distortion interference of the T2 relaxation signal caused by the internal magnetic field gradient generated by paramagnetic minerals (such as pyrite, chlorite, and magnetite). This mechanism significantly improves the characterization accuracy and reliability of nanoscale micropores.

[0098] This invention employs a stitching and fusion method based on a pore size-dependent dynamic weighting function. It assigns high weight to SEM pore size distribution in the nanometer range (utilizing SEM's direct geometric resolution) and high weight to CT pore size distribution in the micrometer range (utilizing CT's three-dimensional volumetric statistical advantages). In overlapping regions, it uses the corrected NMR pore size distribution as the fluid volume conservation benchmark for a smooth transition, achieving seamless and continuous characterization from nanometer (<100nm) to micrometer (>1μm) and then to the macroscopic scale. Simultaneously, based on the pore type classification results from image segmentation, it dynamically assigns differentiated surface relaxation rates to organic pores, inorganic mineral pores, and microfractures. It uses cross-attention weights to calculate weighted equivalent surface relaxation rates to perform differentiated pore size conversion, making the pore size mapping more consistent with the physical reality of different pore types in shale. This provides a high-precision quantitative basis for evaluating the oil-bearing potential of shale reservoirs and designing development plans.

[0099] This invention addresses the engineering challenge of the order-of-magnitude difference in physical spatial observation scale between SEM images (typically with fields of view ranging from a few micrometers to tens of micrometers) and Micro-CT images (typically with fields of view ranging from hundreds of micrometers to millimeters). It employs a hierarchical heterogeneous graph structure: a nanoscale pore topology subgraph is constructed based on SEM images, and a micrometer-scale pore topology subgraph is constructed based on Micro-CT images. Rigid physical spatial edge connectivity is not forcibly established between the two subgraphs; instead, statistically significant high-dimensional feature interactions are achieved through cross-scale aggregation operators. The global pooling feature summary of one scale subgraph is injected as additional contextual information into the node feature update process of the other scale subgraph. This design overcomes the engineering bottleneck of achieving absolute alignment of microscopic coordinates during sample preparation and scanning, and effectively fuses cross-scale statistical features while maintaining the topological independence of each scale subgraph.

[0100] The above are merely preferred embodiments of the present invention and are not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for cross-scale porosity characterization of shale reservoirs based on graph neural networks, characterized in that, Includes the following steps: Collaborative acquisition of low-field nuclear magnetic resonance raw spin echo data and high-resolution scanning electron microscopy or micron-CT image data of shale samples from the same stratum; The image data is segmented to create a pore topology map structure. Relaxation features are extracted from the raw nuclear magnetic resonance spin echo data, and an initial inversion is performed to obtain the initial T2 distribution; The pore topology graph structure data, relaxation features, and initial T2 distribution are input into a graph neural network multimodal fusion model. Geometric features and relaxation features are aligned in the shared hidden layer space to generate an adaptive correction function. This includes: establishing a cross-domain mapping between the T2 distribution and pore nodes through a cross-attention mechanism; using the embedding representation q(m) of the T2 components as the query vector; and applying this mapping to the pore node embedding set { output by the graph neural network}. Attention-weighted aggregation is performed to obtain a comprehensive geometric feature representation. ,in, , i Indicates the pore node index. N Indicates the total number of pore nodes; The adaptive correction function is expressed as follows: ； ； in, For feature dimension, This represents a correction network consisting of fully connected layers, z nmr This represents the latent space eigenvectors derived from the NMR relaxation branch. This represents the learnable parameters of the calibration network. f initial ( m The magnitude of the m-th T2 component obtained from the initial inversion is given by [value]. A non-negative truncation operation is applied to the output of the correction network to ensure [value]. ; The signal distortion of the initial T2 distribution in the low signal-to-noise ratio range is nonlinearly reconstructed using the adaptive correction function to obtain the corrected T2 distribution. Based on the corrected T2 distribution, a full-scale pore size distribution model and multiphase fluid saturation are output.

2. The method for cross-scale porosity characterization of shale reservoirs based on graph neural networks according to claim 1, characterized in that: When performing aperture segmentation on the image data, a U-Net network with a fusion attention mechanism is used. The U-Net network includes an encoder path and a decoder path. The encoder path extracts multi-scale features through multi-level downsampling, and the decoder path fuses features of different scales through skip connections. The attention mechanism introduces a spatial attention gating unit in the skip connections. The segmentation output includes a classification label map of organic matter pores, intergranular pores, and microcracks.

3. The method for cross-scale porosity characterization of shale reservoirs based on graph neural networks according to claim 1, characterized in that, The data used to construct the pore topology map includes: Each segmented independent pore body is abstracted as a node, and the node feature vector includes at least specific surface area, pore volume, aspect ratio, perimeter and equivalent diameter; the pore throats connecting adjacent pore bodies are abstracted as edges, and the edge feature vector includes at least pore throat radius and connectivity weight; and an adjacency matrix is ​​constructed to describe the topological connection relationship between pore bodies.

4. The method for cross-scale porosity characterization of shale reservoirs based on graph neural networks according to claim 1, characterized in that, The initial inversion adopts a multi-exponential decay model based on Tikhonov regularization, and the inversion objective function is expressed as: ； The constraints are as follows: A represents the kernel matrix. Let represent the T2 distribution vector to be determined, M represent the measurement echo vector, λ be the regularization parameter, and L be the regularization operator.

5. The method for cross-scale porosity characterization of shale reservoirs based on graph neural networks according to claim 1, characterized in that: The extraction of the relaxation features is achieved through a one-dimensional convolutional neural network, which consists of three stages: convolution operation, ReLU nonlinear activation, and max pooling downsampling. After being cascaded through multiple layers, the compressed relaxation feature vector is output.

6. The method for cross-scale porosity characterization of shale reservoirs based on graph neural networks according to claim 1, characterized in that, The graph neural network multimodal fusion model adopts a graph attention network architecture, and the attention coefficient is calculated using the following formula: ； The node feature update rule is as follows: , , Representing node characteristics, Representing edge features, W, W e Represents the learnable weight matrix. Let N(i) represent the attention parameter vector, and N(i) represent the set of neighbor nodes. This represents a non-linear activation function.

7. The method for cross-scale porosity characterization of shale reservoirs based on graph neural networks according to claim 1, characterized in that: In the shared hidden layer space, feature alignment is performed through a contrastive learning loss function, so that the geometric features and relaxation features from the same layer and similar lithofacies are close in the hidden space, while features from different lithofacies or with significantly different pore structures are far away in the hidden space. The contrastive learning loss function adopts the InfoNCE loss, expressed as: ； Among them, z geo Let represent the latent space feature vector from the topological branch of the image, and sim(·,·) denote the cosine similarity function. This indicates temperature hyperparameters. V This indicates the number of negative samples compared.

8. The method for cross-scale porosity characterization of shale reservoirs based on graph neural networks according to claim 1, characterized in that, The full-scale pore size distribution model adopts a splicing and fusion method based on a pore size-dependent dynamic weight function, which is expressed as follows: ； Among them, w nano (r), w micro (r) and w nmr (r) are all adaptive fusion weighting functions that dynamically change with the pore radius r, and the weighting functions satisfy... + + =1, PSD SEM (r) represents the nanoscale pore size distribution based on SEM images; PSD CT (r) represents the micrometer-scale pore size distribution based on CT images; This serves as a benchmark for the full-scale continuous pore size distribution obtained by differential surface relaxation rate conversion based on the corrected T2 distribution; Based on the pore type classification results from image segmentation, differentiated surface relaxation rates are assigned to organic pores, inorganic mineral pores, and microcracks. Utilizing cross attention weights β m,i Calculate the weighted equivalent surface relaxation rate ,implement Aperture transformation mapping; In the formula, , k Indicates the pore type, p k (m) This represents the expected proportion of pore type k in the m-th T2 component.

9. A graph neural network-based multi-scale porosity characterization system for shale reservoirs, used to implement the graph neural network-based multi-scale porosity characterization method for shale reservoirs as described in any one of claims 1 to 8, characterized in that, include: The data acquisition module is used to collaboratively acquire low-field nuclear magnetic resonance raw spin echo data and high-resolution scanning electron microscope or micron CT image data of target shale samples in the same stratum; The image segmentation module uses a U-Net network with an attention fusion mechanism to perform pore space segmentation processing on high-resolution image data, and identifies and classifies organic matter pores, intergranular pores and microcracks; The topology graph construction module constructs porous topology graph structure data containing node feature matrices and edge feature matrices based on the segmentation results; The NMR feature extraction module is used to extract features from nuclear magnetic resonance echo data to obtain relaxed feature vectors, and to perform initial T2 inversion based on Tikhonov regularization to obtain the initial T2 relaxation distribution sequence. The GNN fusion correction module aligns geometric and relaxation features in the shared hidden space using a graph neural network multimodal fusion model, generating an adaptive correction function to nonlinearly reconstruct and correct signal distortion in the low SNR region based on the initial T2 relaxation distribution. Generating the adaptive correction function involves: establishing a cross-domain mapping between the T2 distribution and the aperture nodes through a cross-attention mechanism; using the embedding representation q(m) of the T2 components as the query vector; and applying the aperture node embedding set {...} output by the graph neural network. Attention-weighted aggregation is performed to obtain a comprehensive geometric feature representation. ,in, , i Indicates the pore node index. N Indicates the total number of pore nodes; The adaptive correction function is expressed as follows: ； ； in, For feature dimension, This represents a correction network consisting of fully connected layers, z nmr This represents the latent space eigenvectors derived from the NMR relaxation branch. This represents the learnable parameters of the calibration network. f initial ( m The magnitude of the m-th T2 component obtained from the initial inversion is given by [value]. A non-negative truncation operation is applied to the output of the correction network to ensure [value]. ; The results output module is used to output the full-scale pore size distribution model and the quantitative evaluation results of multiphase fluid saturation.

Images