Seismic data is a type of sequence data, and the attention mechanism is inspired by the selective attention mechanism of human vision, which focuses on critical information while ignoring secondary information. It allows sifting through complex data to find helpful information for the present. The core idea of the attention mechanism is to selectively focus on input information by assigning weights that filter out less important information from a large amount of data. The multi-headed attention mechanism is the key component of the Transformer architecture, entirely based on the attention mechanism. It uses multiple copies of the single-headed attention mechanism to extract different information, and the outputs of these heads are concatenated and passed through a fully connected layer to produce the final output.

The multi-headed attention mechanism is defined as follows (Vaswani et al., 2017): Q , K , V = P Q K V z t (1) A t t e n t i o n Q , K , V = s o f t m a x Q K T D V (2)

Where Q denotes the information query set, K is the corresponding similarity set, and V denotes the value set. The matrices of Q , K and V are projected by the inputs in the attention module using the fully connected layer shown in Eq 1. The scaling factor 1 / D is applied to the parametric gradient, the softmax function such that the sum of Q K T / D is equal to 1. The attention weight is denoted as W = s o f t m a x Q K T / D , where W is based on the pairwise correlation between the query set Q and the key set K , D denotes the dimension of the sequence. In addition, the attention weights allow us to extract relevant features using the value set V. For more details on the Transformer architecture, refer to Vaswani et al. (2017).

When using the LSTM network alone to input a long sequence for prediction, the gradient update may decay quickly, hindering the update of sequence data to some extent and making it challenging to represent the feature vector effectively. Additionally, subsequent input data also overwrites the previous input, resulting in the loss of details. However, by using the multi-head attention module and performing a linear transformation to generate Q , K , and V matrices, after expanding the original low-dimensional information to a higher dimension, taking into account the differences between seismic data and natural language, and at the same time having the interaction capability of global information for the input of long sequence multiple seismic attribute data. It not only enhances the influence of strong signals of seismic attribute data but also focuses on weak signal anomalies, avoiding the problem of losing important feature information and reducing the interference of unnecessary information in the sequence. The joint LSTM can effectively learn the rich features and laws in the input data to obtain the porosity data in dense reservoirs. Therefore, the attention mechanism is auxiliary in extracting complex data features.

According to the Temporal Fusion Transformer architecture (Lim B et al., 2021), some modules were modified accordingly for the large seismic data volume. We added a gating and variable selection network to the Transformer for seismic data.

In geological studies, lithology typically varies with depth Logging data provides information about the formation rocks, while seismic attributes reflect different features of the formation. The corresponding values for porosity prediction should be a weighted sum of adjacent characteristic responses that have a certain correlation with each other. Therefore, when we establish the relationship between porosity and external input of various seismic attributes, we should consider the local correlation, the trend of change with depth, and the relationship of adjacent strata. However, since the relationships between external seismic attribute inputs are unknown and uninterpretable for porosity prediction, it is difficult to determine which variables are correlated and the degree of nonlinear processing required. In summary, we need a module for the nonlinear processing of the inputs in this network, so we added a gated residual network (GRN) as a building block for our network. The GRN would accept a primary input a and context vector c , giving it the functionality to obtain more accurate and comprehensive features. G R N ω a , c = L a y e r N o r m a + G L U ω η 1 (3) η 1 = W 1 , ω η 2 + b 1 , ω (4) η 2 = E L U W 2 , ω a + W 3 , ω c + b 2 , ω (5)

Where ELU denotes the unit activation function (Clevert, Unterthiner, & Hochreiter, 2016), η 1 ∈ R d m o d e l , η 2 ∈ R d m o d e l are intermediate layers, L a y e r N o r m is a standard layer normalization by Lei Ba, Kiros and Hinton (2016), ω is the metric indicating weight sharing, b · ∈ R d m o d e l are the weights and biases. The ELU activation function will then function as a recognition function when W 2 , ω a + W 3 , ω c + b 2 , ω is greater than zero, and ELU activation will produce a constant output when W 2 , ω a + W 3 , ω c + b 2 , ω is less than zero. We use a component gating layer based on Gated Linear Units (GLUs) (Dauphin et al., 2017) to provide flexibility to suppress certain repetitions of unwanted parts of a particular seismic attribute dataset.

The variable selection also allows TP to remove any unnecessary noisy inputs that may negatively affect the prediction. However, identifying null or invalid values of the seismic attribute data manually in a large work area can be time-consuming and laborious. Hence, we introduced a variable selection module (Gal & Ghahramani, 2016) to automate this process. A linear transformation of the variables is also applied to convert each post-input variable into a d-dimensional vector, which matches the dimensions in the subsequent layers and is mainly used to achieve jump connections. The weights vary for different data (determined by the correlation of each attribute with the porosity). The variable selection network on the seismic data input is presented below, as shown in Figure 2. Note that the variable selection network is the same for the other data set inputs.

We denote X i j ∈ R d m o d e l as the j t h variable input value in a time-domain earthquake sequence, where Ξ i = X i 1 T , … , X i m T T is the flattened vector of which all have been entered. The weights for variable selection are then generated by feeding both Ξ i and additional context vectors c into GRN and then into the Softmax layer: υ χ t = S o f t m a x G R N υ χ Ξ t , c (6)

Where υ χ t ∈ R m χ is the vector of variable selection weights and c is obtained from the static covariate encoder (Section 2.3).

An additional non-linear processing layer is applied during the input process by passing each sequence element through its own GNR. X ∼ i j = G R N X ∼ j X i j (7) where X ∼ i j is the feature vector after processing the seismic attribute j . We note that each sequence has its own G R N X ∼ j , and the weights are shared globally. The processed features are then weighted and combined according to their variable selection weights as follows: X ∼ i = ∑ j = 1 m υ χ t j X ∼ i j (8)

Where υ χ t j is the j t h element of vector υ χ t .

In summary, the correlation between seismic attribute data and porosity is unknown, and adding the GRN module provides the advantage of complementing the seismic attribute preference. Additionally, because the input of various seismic attribute data may contain many invalid values, the variable selection network is used to further optimize the input data and prevent significant errors in the prediction results.

Due to the variation of rock properties with depth and the correlation between seismic response adjacencies, it is crucial to consider not only the local correlation between seismic attribute data but also the trend of seismic data with depth and adjacency information when establishing the mapping relationship between different attribute data and porosity. The TP model can integrate information from the logging porosity data. A separate GNR encoder is then used to generate a context vector. The context vectors are connected to different locations in the network architecture where static variables play an essential role. These include the selection of data for additional seismic attributes, the underlying processing of features, the use of logging curves to guide the identification of the context vector, and the enrichment of features for data adjacency.

We have considered the use of the LSTM encoder-decoder as a building block for our prediction architecture, following the success of this approach in typical sequential coding problems (Wen et al., 2017; Fan et al., 2019; Wang et al., 2022). Due to the specificity of seismic data, the critical points in the data we apply are usually determined based on the surrounding values. Therefore, we propose to use a sequence-to-sequence layer to process the data naturally. Send X i − k : i to the encoder and X i + k : i + τ max to the decoder. A set of uniform features is generated as input to the decoder itself, using ϕ i , n ∈ ϕ i , − k , … , ϕ i , τ max , and n is a position index. In addition, to allow static metadata to affect local processing, we use a context vector from a static covariate encoder to initialize the cell state and hidden state of the first LSTM in the layer, respectively. We also employ a gated jump connection on this layer at: ϕ ∼ i , n = L a y e r N o r m X i + n + G L U ϕ ∼ ϕ i , n (9)

Where n ∈ − k , τ max is a location index, i denotes the i t h sequence data.

Since static covariates often have an enormous impact on features, we introduce a layer simultaneously to augment the corresponding features with static metadata. For a given location index n , the form of: θ i , n = G R N θ ϕ ∼ i , n , c (10)

Where the whole layer shares the weights of G R N θ , and c is the context vector from the static covariate encoder. Finally, the attention layer, the decoder mask (Vaswani et al., 2017; Li et al., 2019b), is applied to the multi-head attention mask layer to ensure that each dimension notices its previous features.

Considering the completeness and universality, many seismic attributes need to be extracted, resulting in redundant information. Seismic attribute optimization is to preferably select the attributes with better correlation with the target reservoir parameters from a large number of attributes to eliminate the redundant information and thus improve the prediction accuracy. The preferred effect of seismic attributes is often reflected in oil and gas prediction results. In this study, the analysis of the rendezvous plots was first conducted using the Pearson correlation coefficient (PCC) as the criterion. The PCC reflects the degree of similarity of each unit change between two variables, and the closer the result is to 1 indicates the more significant correlation between the variables, which is calculated by the following equation: P C C = C o v y , y ^ σ y σ y ^ (11)

Where C o v y , y ^ denotes the covariance of the two variables y and y ^ , used to describe the linear relationship between the two variables, and σ is the standard deviation.

This study analyzed porosity separately with multiple attributes, and we used the attribute preference method to perform preliminary screening of each seismic attribute. The multiple seismic attributes shown in Table 1 include post-stack seismic attributes and various derived post-stack attributes in order to avoid the multi-solution nature of single-parameter inversion to initially screen out the useless data with enormous influence. Figure 1 shows a representative rendezvous analysis of several attributes, and other rendezvous analyses are also performed in the same way we can select the corresponding seismic attributes. Finally, 12 seismic attribute bodies such as Q Factor, Fluid Factor, PG, P-impedance, and variance are selected as input according to PCC, as shown in Table 2 and Figure 3.

This study uses three primary metrics to evaluate the performance of porosity prediction, namely, mean square error (MSE), Pearson correlation coefficient (PCC), and coefficient of determination ( R 2 ), to quantitatively evaluate the inversion results.

Mean squared error. MSE is the sum of squares of the corresponding point errors between the predicted and actual data. The smaller the value, the better the fit between predicted and actual data. It is defined as shown in Eq 12. M S E = 1 N ∑ i = 1 n y i − y ⌢ i 2 (12)

Where y i and y ⌢ i are the predicted and actual values of the series, and N is the number of sampling points.

Pearson Correlation Coefficient. PCC indicates the overall fit between the predicted and actual data. Its range is [-1,1]. The larger the value, the stronger the linear correlation. It is defined in Eq 11.

Coefficient of determination. R 2 is used to assess the fit between the variables and considers the mean square error between the predicted and actual data. It ranges from [0,1]. The larger the value, the better the fit between the variables. The following equation gives it: R 2 = 1 − ∑ i = 1 n y i − y ⌢ i 2 ∑ i = 1 n y i − y ¯ i 2 (13)

Where y ¯ is the average value of y i i = 1,2 , . . . , n .

In order to assess the validity of our proposed method, we conducted a method analysis experiment using logging data from two wells, X and Y, in a dense sandstone reservoir in an actual working area. The logging depth of well X was between 2750 m and 2940m, while that of well Y was between 2714m and 2870 m. The logging curves included porosity, density, spontaneous potential, p-wave velocity, natural gamma ray, deep lateral resistivity, compensated neutron, and others. After removing any outliers, we obtained two sets of logging data from wells X and Y for training and test datasets in the Transformer architecture. First, we screened the sensitive parameters of well X using Pearson correlation coefficients. Next, we used the Transformer model built with the logging data of well X to predict the porosity of well Y. To compare our model with other existing models, we also used bi-directional Long Short-Term Memory (BLSTM) and convolutional neural network (CNN) in this experiment.

The correlation between the porosity and other logging parameters was obtained by correlating the petrophysical logging parameters of the two test study wells, as shown in Figure 4.

The logging parameters selected for this study are density (DEN), spontaneous potential (SP), p-wave velocity (VP), natural gamma (GR), acoustic log(DT), and compensate neutron log(CNL). From the correlation graph shown in Figure 4, it is evident that CNL, DEN, DT, and GR are highly correlated with porosity (POR), whereas natural potential and longitudinal velocity exhibit weaker correlations. Therefore, we have used the four input datasets with higher correlation for predicting porosity in tight reservoir logs. In this study, we have used Root Mean Square Error (RMSE), Mean Absolute Error (MAE), and Mean Absolute Percentage Error (MAPE) as evaluation metrics to assess the performance of the model. These metrics measure the deviation between the predicted and actual data, with lower values indicating better performance.

We trained the model on the preferred well X data and used it to predict the porosity of well Y. The results of the prediction are presented in Figure 5. Using deep learning techniques and sensitive feature parameters, such as CNL, DEN, DT, and GR, we obtained good results in porosity prediction of wells. Three models, namely, Transformer, BLSTM, and CNN, were used for porosity prediction, and while errors existed, the overall fit to the actual porosity curve was satisfactory. As illustrated in Figure 5, the predicted values of the Transformer model were closer to the actual values than those of the BLSTM and CNN models, especially in the areas of significant variation, where the results were more satisfactory. The prediction errors are summarized in Table 3, and it can be concluded that the Transformer model has more advantages.

In this study, we applied the model to the actual data of an exploration area in eastern Sichuan for tight sandstone oil and gas exploration. It is difficult to describe the reservoir in this exploration area; the seismic response of the submerged river sand body is not obvious and variable, the reservoir thickness of the river sand body is large, but the porosity is small, the porosity prediction is difficult, and it is a dense gas reservoir. We used the logging data of 5 wells and optimized 12 types of seismic volume attribute data in this work area. We divided our dataset into three parts - a training set for learning, a validation set for hyperparameter tuning, and a test set reserved for performance evaluation. For hyperparameter optimization, we use the random search method. Table 4 shows the entire search range of all hyperparameters as follows and the optimal model parameters.

● State Size–10,20,40,80,160,240

During the training process of our dataset, a significant amount of seismic data input is required. However, using the basic architecture of the Transformer model, the entire architecture can be deployed on a single GPU without consuming many computer resources. For instance, when we used NVIDIA Quadro GP100 GPU for the experimental data, our optimal parameter TP model only required less than 2 hours to complete training. Figure 6 presents the comparison of well porosity prediction results, in which the red curve is the actual value, the blue curve is the CNN prediction value, the purple curve is the BLSTM prediction value and the green curve is the TP prediction value. With the CNN and BLSTM data were added to compare network effects in this study. From the comparison in the figure, we can find that the prediction result of TP is more accurate than CNN or BLSTM, and the mean square error (MSE), Pearson correlation coefficient (PCC), and are used as regression evaluation indicators. Table 5 selects four wells (coded by A, B, C, and D), and it can also be proved from Table 5 that our TP can obtain more accurate results.

Figure 7 shows the porosity result profile predicted by the TP model in this study. Specifically, Figure 7A presents the porosity result profile utilizing standard software for frequency division inversion, whereas Figure 7B represents the porosity inversion of TP training. In the porosity result profile, the red area represents the distribution range of high porosity. The porosity curve of critical well location A in this work area is utilized to verify the prediction results. It is known that the porosity of tight reservoirs inverted by traditional methods is generally small. The important basis is the well-logging porosity, which is significantly affected by the logging data, so the overall resolution is low. Figure 7A also proves that the overall porosity resolution inversion of the tight reservoir is not high, and the lateral continuity is poor. During porosity comparison with logging, the overall porosity in Figure 7A is relatively small due to the complex characteristics of tight reservoirs. However, it can be found from the logging report that the actual porosity range is 0–8%. Hence, the single inversion method still has limitations and does not match the porosity reservoir characteristics of the actual work area. In contrast, Figure 7B illustrates the resulting porosity profile predicted by the TP model. It is observed that the high values are also concentrated near Well A compared to Figure 7A. The overall porosity range also reaches 0–8%, but Figure 7B has a higher resolution. The comparison around 2058 ms shows that the porosity distribution in Figure 7B is more accurate in the lateral direction. Figure 7A shows that the traditional method cannot accurately invert the pore distribution in the right half, far away from the well location, due to the limitation of logging data. Table 6 shows the well and porosity profile prediction and the MSE, PCC, and R 2 of logging and input porosity profile, respectively.