Permeability extraction from multiple well logs using particle swarm optimization based factor analysis

In this paper, we present an innovative factor analysis algorithm for hydrocarbon exploration to estimate the intrinsic permeability of reservoir rocks from well logs. Unlike conventional evaluation methods that employ a single or a limited number of data types, we process simultaneously all available data to derive the first statistical factor and relate it to permeability by regression analysis. For solving the problem of factor analysis, we introduce an improved particle swarm optimization method, which searches for the global minimum of the distance between the observed and calculated data and gives a quick estimation for the factor scores. The learning factors of the intelligent computational technique such as the cognitive and social constants are specified as hyperparameters and calculated by using simulated annealing algorithm as heuristic hyperparameter estimator. Instead of the arbitrary fixation of these hyperparameters, we refine them in an iterative process to give reliable estimation both for the statistical factors and formation permeability. The estimated learning parameters are consistent with literature recommendations. We demonstrate the feasibility of the proposed well-log analysis method by a Hungarian oilfield study involving open-hole wireline logs and core data. We determine the spatial distribution of permeability both along a borehole and between more wells using the factor analysis approach, which serves as efficient and reliable multivariate statistical tool for advanced formation evaluation and reservoir modeling.


Introduction
Wireline logs including abundant in situ measurement information taken from boreholes can be effectively used to derive basic petrophysical parameters of hydrocarbon reservoirs (Serra 1984). Most of these formation properties (e.g., porosity, shale volume, water saturation etc.) are usually estimated in reliable quality-checked inversion procedures (Ball et al. 1987;Alberty and Hashmy 1984;Dobróka et al. 2016). However, permeability of rock formations cannot be directly estimated by the inversion of conventional well logs, since the available tool response equations do not include it explicitly. Therefore, it is generally estimated empirically by using some types of well logs (Wyllie and Rose 1950;Timur 1968;Coates and Dumanoir 1974) and these estimations often need to be calibrated by core data (Singh 2018). Due to different methods and conditions, the estimation error of permeability is rather large compared to the other petrophysical parameters (e.g., porosity and shale volume) that are required for evaluating the productivity of hydrocarbon reservoirs. Exploratory data analysis methods including regression analysis, dimension reduction and pattern recognition techniques are widely used in geosciences for different purposes. Tang and White (2008) made a comparative study on several statistical tools for facies classification. Jafarpour and McLaughlin (2009) applied a discrete cosine transform for reservoir characterization. Soares et al. (2020) used dictionary learning for the sparse representation of 4D seismic data to highlight only the main features of the dataset. Cosultchi et al. (2012) utilized principal component analysis and cluster analysis for aquifer identification in petroleum reservoirs. Generally, some linear regression analysis approach is used for relating the permeability to existing petrophysical parameters (Dahraj and Bhutto 2014;Handhal 2016). Multivariate statistical tools have been applied for permeability estimation in Al-Mudhafar (2020).
We present a new alternative to predict permeability of hydrocarbon reservoirs, which is based on a metaheuristic approach for factor analysis. The simultaneous processing of well logs allows not just to reduce the dimensionality of the data space, but to find correlation between the extracted factor variables and reservoir parameters. As one of the most important parameters, we estimate the intrinsic permeability utilizing all available well logs leading to more reliable estimates. Szabó and Dobróka (2018) suggested a statistical approach for the processing of well logs, in which a differential evolution algorithm based factor analysis was introduced for shale volume estimation. For solving the problem of factor analysis, we currently apply the particle swarm optimization (PSO) method proposed by Kennedy and Eberhart (1995), which is less computationally expensive and requires less parametrization than genetic algorithm (Rahmat-Samii 2003;Donelli et al. 2006). This metaheuristic technique is based on swarm intelligence (Del Ser 2019), which is widely used in engineering applications. For instance, Nguyen et al. (2020) utilized it for feature selection in data mining. Self et al. (2016) suggested PSO for selecting optimal drilling parameters to minimize the overall cost. Shaw and Srivastava (2007) evaluated the applicability of PSO to inversion of direct current, induced polarization and magnetotelluric data and concluded that the results are consistent with ridge regression and genetic algorithm based inversions. Essa and Elhussein (2018) presented a geophysical application, in which magnetic anomalies were interpreted successfully with the help of this metaheuristic.
Factor analysis is a multivariate statistical tool that is generally used to decrease the dimensionality of multi-parameter problems, while keeping most of the information contained in the original variables to simplify the data interpretation. The newly derived factors are necessarily uncorrelated and may reveal hidden information from the measurements (Lawley and Maxwell 1962). Thus, for processing wireline logging data it can be effectively applied because the measured datasets are usually quite large due to several applied logging tools and numerous measured depth points. The extracted factors from wireline logs are referred to as factor logs that can be linked to petrophysical parameters through regression analysis (Szabó 2011). In this way, factor analysis processes all available well logs to give an estimation to these parameters, while for instance, simple deterministic methods normally rely only on one log type to derive one parameter at a time. Since it is an effective approximate method, it can also process data for which an exact mathematical relation with petrophysical properties is not available (e.g., the caliper log), however they may contain useful lithological or physical information.
We show a two-level embedded computing algorithm for improving the performance of factor analysis. In the inner loop, the PSO algorithm runs to estimate the factor scores while the factor loadings are initially derived from the non-iterative approach of Jöreskog (2007). In the outer loop, the hyperparameters of the same PSO procedure called the cognitive and social parameters are refined by a Metropolis criterion based simulated annealing method as global optimization tool avoiding the possible localities of the objective function (Metropolis et al. 1953). These control parameters are important learning factors, which are normally chosen as arbitrary constants at the beginning of the PSO procedure. In this study, we suggest a much more intelligent manner as estimating them automatically. This approach was first proposed for the interpretation of direct push logs to estimate the volumetric soil moisture content of the shallow unsaturated zone (Abordán and Szabó 2021) as a further improved factor analysis technique used in Szabó et al. (2012). In this study, the factor scores are crossplotted with reservoir properties. We assume the extracted factors hold quantitative information on the permeability of deep reservoir rocks. We extend the factor analysis algorithm to process multiwell logging data. Our permeability prediction results are shown involving single and multi-borehole data collected from a Hungarian gas field.

Traditional method of factor analysis
To initialize factor analysis, we build a common data matrix (D) from the standardized data of S number of measured well logs, in which every column represents a dataset recorded by different logging tools, while N number of rows contain all measured depths along the processed interval Factor analysis is based on the following model where F is the factor scores matrix of size N-by-R (R is the number of derived factors). Quantity L is the factor loadings matrix of size S-by-R and E is the residual error matrix of size N-by-S (symbol T is the matrix transpose). According to Eq. (2), the observed wireline logs are expressed as a linear combination of the factors. The factor loadings in matrix L show the partial correlations between the observed and derived statistical variables. The largest part of total variance represented by matrix D is explained by the first factor, the scores of which are found in the first column of matrix F. Given that the factors are uncorrelated, we write the correlation matrix of standardized data as where denotes the diagonal matrix of error variances independent of those of the common factors in matrix F. The calculation of correlation matrix R is easy to handle if matrix E is neglected in Eq.
(2). This approximation helps to avoid the estimation of specific variances by cancelling matrix in Eq. (3). Since the part of the total variance that may not be described by the common factors is neglected, this approach leads to an approximate solution for the factors using the popular method of Principal Component Analysis (PCA). Although PCA does not account for residual matrix E, but it gives a quick and unique (algebraic) solution for extracting the principal components (Maćkiewicz and Ratajczak 1993).
We approximate the factor loadings non-iteratively by a scaling-free estimation method (Jöreskog 2007) where denotes the matrix of the first R number of eigenvalues of the empirical covariance matrix , and the corresponding R number of eigenvectors are given in matrix , and U is a freely selected R-by-R orthogonal matrix. Quantity I denotes the identity matrix and θ is a properly chosen constant that must be accepted smaller than 1. In traditional dimension reduction approaches, the factor scores are calculated by the maximum likelihood method. By assuming the hypothesis of linearity gives an unbiased estimate to the factor scores (Bartlett 1937). Since Eq. (5) is analogous to a least squares solution weighted by the specific variances (Menke 2012), the localities of the objective function to be minimized cannot be effectively avoided by the gradient-based searching for the absolute minimum.

Particle swarm optimization based factor analysis
To improve the performance of data reduction, we suggest PSO as a global optimization procedure to determine the vector of factor scores. Abordán and Szabó (2018) established the basic principles of the PSO assisted factor analysis (FA-PSO). To give a metaheuristic estimation of the factors, we reformulate the model of factor analysis given in Eq.
(2) as where d (m) is the vector of standardized wireline logging data of size SN,L denotes the matrix of factor loadings of size NS-by-NR, f is the expanded vector of factor scores of size RN and e is the vector of approximation errors of size SN (Szabó and Dobróka 2017). As a first step, we determine matrixL by Eq. (4) and rotate it using the varimax technique for easier interpretation (Kaiser 1958). We give the optimal values of the factor scores f by using PSO, which can be considered as an iterative inversion procedure. We use the following L 2 -norm based function as objective function where d (c) contains the standardized calculated wireline logging data. Here, calculated data refers to the termLf from Eq. (6), which also gives the possibility for calculating theoretical well logs that can be regarded as the solution of the forward problem of the inversion procedure. Previous researches suggest that it can be even used to generate theoretical well logs for intervals where measurements are missing, due to the high correlation of certain log types to the first factor (Szabó et al. 2012). By the minimization of the objective function in Eq. (7), we update the factor scores by PSO, while we keep the factor loadings constant to reduce runtime. It must be mentioned that it is possible to simultaneously refine the factor scores and loadings using the well-known maximum likelihood method, or an iteratively reweighted factor analysis for seeking a more robust solution (Szabó and Dobróka 2017). Particle swarm optimization originally introduced by Kennedy and Eberhart (1995) was developed to mimic the social behavior of animals. Generally, it was applied for its easy implementation and high efficiency in solving complex problems (Kumar et al. 2017). To find the optimal solution, PSO applies a swarm of particles generated randomly by following uniform distribution in a pre-defined search space. Given an GEM -International Journal on Geomathematics (2022) 13:10 N-dimensional search space, vector f i = (f i1 ,f i2 ,…,f iN ) T denotes the position of the ith particle. In our problem, quantity f i denotes the ith possible solution for the vector of factor scores. We store the particle's velocity in vector v i = (v i1 ,v i2 ,…,v iN ) T , which defines its movement in each iteration within the search domain. The swarm of particles seek the minimum of the objective function defined in Eq. (7) by utilizing the following equations where L denotes the number of particles in the swarm (i = 1,2,…,L), t stands for iterations (t = 1,…,t max ) where t max is the final iteration step. During the search for the global optimum, PSO utilizes a memory which is updated after each iteration. Every particle's best position is held in p i = (p i1 ,p i2 ,…,p iN ) T , while vector g stores the best position found by the whole population. Kennedy and Eberhart (1995) suggested to set the cognitive and social constants c 1 and c 2 equal to 2 as default, separately, during the initialization. These control parameters define to what extent the personal best position of particles p i and the very best position of the swarm g affect the movement of the ith particle during the search. In Eq. (9), r 1 and r 2 are uniformly distributed random parameters from 0 to 1, while w represents an inertia weight that was suggested by Shi and Eberhart (1998) to more effectively control the optimization procedure.
To start the FA-PSO procedure for finding the optimal vector of factor scores f, we generate a population of particles with uniform distribution within the search space, which can be pre-defined by solving Eq. (5) for the factor scores. As for w, we use a chaotic descending inertia weight scheme to better control the velocity of particles (Feng et al. 2007). We perform it in three steps: first we generate a random number (z) from the range 0 to 1. Then we set z by logistic mapping according to z = 4z(1-z) and finally where w 1 and w 2 represent the starting and the final values of inertia weight during the iterative procedure, respectively. Then in each iteration step, we re-determine the particles' positions according to Eqs. (8-9), and re-calculate the objective function defined in Eq. (7) with the new factor scores f. We continuously update the particles until the maximum number of iterations is reached, thereafter we accept the particle with the lowest data misfit as the solution for the factor scores.

Permeability estimation by factor analysis
Previous studies show that the first statistical factor correlates well with the shale content of various geological formations (Szabó 2011;Szabó and Dobróka 2013;Puskarczyk et al. 2015;Puskarczyk 2020;Abordán and Szabó 2020). Advanced mathematical approaches have been developed to extract the shale volume from well logs as efficient as possible. The iteratively reweighted factor analysis gives a reliable estimation to the factors independent of the statistical distribution of input data (Szabó and Dobróka 2017). In addition to give a robust solution, the calculation of factor scores was established on global optimization by using evolutionary computation (Szabó and Dobróka 2018). By analyzing multiple well logs, we take all the measured information into account during the interpretation, which makes a step forward in increasing the reliability of estimation compared to single depth-point approaches using one or few well logs. Under taking all measured information we mean that we select all suitable well logs as input being highly sensitive to shale volume, in other words only those data types are processed that are significantly influenced by the variation of the relevant petrophysical parameter along the borehole. Based on the above considerations, the first statistical factor extracted from multiple well logs can be interpreted as a lithology indicator, which holds valuable information on the shale content of sedimentary formations. When applying factor analysis to well logs, a strong exponential connection between the first factor (F 1 ) and shale volume (V sh ) can be set in geological formations of different ages. In the Hungarian and American examples given in Fig. 1, the consistence of the results is also confirmed by the validity of the same regression function (with approximately the same regression coefficients and by using proper scaling to the factor scores).
We assume a strong correlation between the absolute permeability and the first principal factor in sedimentary formations. This is partly due to the fact that permeability is to some extent inversely proportional to shale content of clastic formations (Schön and Georgi 2003;Revil and Cathles 1999). In our test site, an empirical relation is found between the decimal logarithm of permeability and shale volume both derived from wireline logs as well as core information (Fig. 2). The correlation coefficient for the first case is − 0.91, while for the latter it is − 0.72, both showing a clear relationship between the two quantities. In well log analysis, for estimating the absolute permeability (K) from wireline logging datasets, usually deterministic approaches are used which rely on porosity (Φ) and irreducible water saturation (S w,irr ) data. The most frequently applied manner is based on the use of similar empirical approaches like the formula of Timur (1968) Analogous studies show that the statistical factors include not only the information on lithological characteristics but also the hydraulic properties of reservoir rocks. For instance, Szabó (2015) explored the hydraulic conductivity of fractured groundwater formations successfully by factor analysis of hydrogeophysical logs. In this study, we test the proposed factor analysis approach to estimate permeability in six wells (Well 1-6) drilled in the Pannonian basin in East Hungary. We investigate unconsolidated upper Miocene gas-bearing formations with high average porosity showing a variability of different grain sizes and lithology as marl, silt and sands. Accordingly, we show that the well log of the first factor derived by the improved FA-PSO approach scaled into the range 0 to 1 can be directly used to estimate permeability along arbitrary GEM -International Journal on Geomathematics (2022) 13:10

Fig. 1
Regression relationships between the first factor and shale volume derived by well log analysis for different hydrocarbon fields intervals. For this purpose, we use a non-linear relationship, which connects the first factor to the intrinsic permeability of hydrocarbon reservoirs. In the given area, the first (properly scaled) factor extracted from the wireline logging dataset by FA-PSO is connected to the decimal logarithm of permeability (K) in the form as where a, b, c, d are area specific coefficients estimated by regression analysis using for instance core data and other preliminary information. The above non-linear relationship is applied to different wireline logging datasets in the investigated hydrocarbon field to demonstrate the feasibility of the proposed method.
Our FA-PSO method is a fully metaheuristic approach, therefore its output is highly affected by the control parameters (i.e., learning factors) that are worth to be set as optimal as possible at initialization. The optimal parameter selection of PSO has been studied by several researches (Jiang et al. 2007;Fernández Martínez et al. 2010;Fig. 2 Empirical connections between the intrinsic permeability and shale volume of Miocene hydrocarbon formations in the investigated site, Pannonian basin, East Hungary Harrison et al. 2019). However, since the optimal control parameters highly depend on the optimization task, here an attempt is made to automatically select the optimal values of parameters c 1 and c 2 of PSO while solving factor analysis. These parameters are usually chosen according to empirical suggestions, which can make a problem ambiguous. To offer a more advanced solution of the presented PSO-based factor analysis, we propose a further developed algorithm to improve the FA-PSO method in regard of the c 1 and c 2 hyperparameters. For this purpose, we search for these parameters by Simulated Annealing method (SA) suggested by Metropolis et al. (1953) in an outer loop of our algorithm in an automated way. Then, in the inner loop, FA-PSO is run with the newly estimated c 1 and c 2 parameters. According to the previously outlined workflow, the procedure starts with the determination of the factor loadings according to Eq. (4). Then, we optimize the factor scores with the help of PSO by minimizing the objective function defined in Eq. (7). However, the control parameters c 1 and c 2 of PSO are also optimized during the procedure by SA in an outer loop as shown in Fig. 3. In the regression part of the workflow, we can correlate the first factor to available permeability values (and other petrophysical parameters). In the prediction phase, we estimate the permeability to missing logging intervals directly from the factor scores using the regression model in Eq. (12).

Single borehole application
We test the feasibility of factor analysis using in situ well logging data in a hydrocarbon well (Well-6). For permeability estimation, we simultaneously process the caliper (CAL), natural gamma-ray intensity (GR), potassium-thorium product (K × TH), spontaneous potential (SP), deep resistivity (R d ), bulk density (ρ b ), neutron-porosity  Table 1. The first factor is strongly related to the natural gamma-ray intensity, potassium-thorium index and bulk density logs that make it a good lithological indicator, which is consistent with previous applications (Szabó and Dobróka 2017). The second factor log correlates with the neutron-porosity and acoustic logs, while the third factor is greatly sensitive to the variation of the deep resistivity log.
Once the factor loadings are determined, we optimize the factor scores by the algorithm of PSO. First, we give an initial solution for the factor scores by solving Eq. (5), thus we define the search space of the PSO process as [−12, 12]. We estimate the amount of variance explained by the factors based on the singular value decomposition of the reduced covariance matrix of standardized data R-= LL T . Since three factors explain more than 90% of total variance of the input data, we extract three factors from the dataset. Thus, the number of unknowns is 4,353 (3 factor × 1451 depth points). For estimating their optimal values by PSO, we initialize 90 particles  (9) to minimize the objective function in Eq. (7) in 5,000 iteration steps. In our experiment, we treat the cognitive and social parameters of PSO as unknown. Our algorithm allows changing them freely, instead of selecting them as arbitrary constants. We set the initial value of c 1 and c 2 equal to 1, respectively, while parameter w is set according to Eq. (10). At the end of the process, the estimated values of hyperparameters are c 1 = 1.44 and c 2 = 2.06. The variation of the control parameters of PSO and data distance as functions of iteration steps is plotted in Fig. 4. At the end of the optimization phase, the data distance calculated by Eq. (7) reaches 0.53. Both plots show a convergent and stable optimization procedure. The independence of results on the random initial particle positions in PSO algorithm has been earlier investigated in the context of inverse problem solutions (Abordán and Szabó 2020). To find the regression relationship between the first factor found by the FA-PSO method and decimal logarithm of absolute permeability, we properly scale the former into the range of 0 to 1. For regression analysis, we use permeability data deterministically derived from Eq. (11). The regression coefficients in Eq. (12) are found to be a = 5.90, b = 1.27, d = 2.28 while c is fixed at 1, because it shows no sensitivity in this case. The correlation relationship between the scaled first factor and absolute permeability is plotted in Fig. 5. The Spearman's rank correlation coefficient of -0.77 between the first factor and the decimal logarithm of permeability indicates a significant (non-linear) inverse relationship. The input well logs of the FA-PSO procedure are plotted in Fig. 6 in the first 8 tracks (note that the depth scale represents relative depth coordinates). The last track indicates that the fit between the factor analysis derived permeability (lg(K)_FA-PSO) and that of derived deterministically by Eq. (11) (lg(K)_DET) is consistent along the investigated interval. The CPU time of the FA-PSO procedure for automatically determined hyperparameters is~7 h (~3000 depth points), while for fixed hyperparameters it is 86 s using an octa-core processor workstation.

Test of hyperparameters
The learning parameters of classical PSO are usually set to default values based on recommendations (Kennedy and Eberhart 1995). We suggest an improved PSO algorithm to improve these parameters, which forms the basis of the workflow presented in Fig. 3. We focus on the estimation of parameters c 1 and c 2 to optimize the PSO process. As input, we select the GR, K × TH, SP, R d , ρ b , Φ N , t logs in a 20 m long interval in Well-3. At the beginning of the SA procedure, both parameters are set as 1 and in every iteration step their value is slightly altered by adding a random perturbation (δ) to both values. We generate parameter δ as random number in each Fig. 6 The result of factor analysis in Well-6, Pannonian Basin, Hungary. The input logs of FA-PSO procedure are plotted in tracks 1-8: borehole caliper (CAL), natural gamma-ray intensity (GR), potassium-thorium product (K×TH), spontaneous potential (SP), deep resistivity (R d ), bulk density (ρ b ), neutron-porosity (Φ N ) and acoustic (P-wave) transit time (Δ t ). The scaled first factor log (F 1 ) is plotted in track 9. In track 10, lg(K)_FA-PSO represents the factor analysis derived logarithmic permeability (originally given in mD), while lg(K)_DET shows the deterministic estimation for the same quantity iteration step from the range −δ max to δ max . We also decrease the maximum value of perturbation (δ max ) iteratively using δ max = δ max × ε, where constant ε is smaller than 1. Then with these new hyperparameters (c 1 and c 2 ), PSO optimizes the factor scores in 200 iterations. We define the energy function of SA as Eq. (7). If the energy difference ( E) in two subsequent SA iterations of the algorithm is negative (i.e., PSO reached a lower data distance), then the new values of c 1 and c 2 are accepted and the procedure is continued. If the energy difference is greater than 0 (i.e., the efficiency of PSO got worse with the new parameters), then the new c 1 and c 2 hyperparameters can be still accepted if a random number generated from 0 to 1 is smaller than acceptance probability P a = exp(− E/T ), where T denotes the actual temperature of the search. This mechanism is the basis of the SA algorithm because it enables the search to get out of local minima of the energy function. To ensure the convergence to the global optimum of the objective function, the temperature must be reduced logarithmically as suggested by Geman and Geman (1984). We change the temperature in the qth iteration as T (new) = T 0 /lg(1 + q), where the initial temperature of the artificial system T 0 is experimentally set to 5 × 10 -6 . We run the SA procedure for 100 iteration steps to optimize the values of c 1 and c 2 as shown in Fig. 7. The newly generated control parameters are tested in each SA iteration step by running the FA-PSO method for 200 iterations for three consecutive times. This is particularly advantageous especially at the early stages, where data distances remain rather large due to the suboptimal values of c 1 and c 2 . The mean of the three consecutive FA-PSO runs for each SA iteration can be seen in Fig. 7 (on the right panel). As the figure indicates, the whole procedure was run 200 times to check its stability. After approximately 40 iteration loops, SA successfully optimizes the values of the learning parameters for the FA-PSO algorithm in all 200 separate program runs and data distance defined by Eq. (7) did not decrease any further. This shows the stability of the SA-based hyperparameter estimation method. In Fig. 8, we illustrate all those value pairs of c 1 and c 2 that were tried by SA during the iterations in 200 repeated program runs. Green dots indicate the ones that did not lead to the global optimum, while the blue ones enabled PSO to converge to the optimal solution of factor scores. According to the results, c 1 has less effect on the convergence and can be selected from a wider range (0.25-2.7), while the optimal range for c 2 is more restricted (1.7-2.6) when solving the problem of factor analysis. The fitted regression line indicates a moderate correlation between c 1 and c 2 . Our result agrees well with those of the literature, for instance Cleghorn and Engelbrecht (2014) derived similar convergent regions for the sum of the two parameters. This test justifies the advantage of using SA to find the optimal combination of hyperparameters. The CPU time of the FA-PSO procedure for automatically estimated hyperparameters is~70 s, while for fixed hyperparameters is 0.27 s using an octa-core processor workstation.

Multi-borehole application
We extend the improved factor analysis to multi-dimensional cases, which allows the simultaneous processing of well logs measured from several neighboring boreholes.  (1) . . .
whereL (h) represents the matrix of factor loadings and vector f (h) contains the estimated factor scores along the hth borehole. There is N h number of measured depth points in the hth wellbore, therefore the total number of data processed by factor analysis is N * = N 1 + N 2 + … + N H . In Eq. (13), we determine the matrix of factor loadings in a similar manner as in the previously presented one-dimensional case by Eq. (4), while we approximate the optimal values of factor scores by PSO. We show a 2D application of FA-PSO to derive the permeability variations in a hydrocarbon zone in the same area of the Pannonian Basin. We gathered the wireline logging data collected from five wells (Well-1, -2, -3, -4, -5) into one dataset. The borehole investigated in Sect. 3.1 (Well-6) falls farther from the line of the other wells, thus we excluded it from our calculation. The GR, ρ b , Φ N , t and shallow resistivity (R s ) logs being completely available all depths serve as the input of 2D factor analysis. By combining the dataset of the five wells, the total number of data is 37,525. The wells are located along a 2,200 m long profile. By estimating two factors, the orthogonally transformed factor loadings are given in Table 2. When calculating the factor scores the learning parameters are simultaneously optimized.
The resultant values of learning parameters are c 1 = 1.59 and c 2 = 2.23 (Fig. 9). Due to the increased number of unknowns, in this case the PSO routine requires 25,000 iteration steps with 90 particles. Based on the calculation of the singular values of the reduced covariance matrix of standardized wireline logs, the first factor explains 77% of the total data variance, while the remaining 23% is explained by the second factor. The first extracted factor correlates well with all processed well logs except the P-wave acoustic interval time log, while the second factor is quite the opposite.
We find a strong correlation between the first scaled factor and the decimal logarithm of permeability in the investigated Miocene formations (Fig. 10). The Spearman's rank correlation coefficient is -0.93. Using the proposed model in Eq. (12), we establish the local regression model, where the available permeability values are obtained from deterministic modeling using Eq. (11). The regression coefficients in this case are a = 7.18, b = 4.5, c = 1.2 and d = -3.8. By interpolating the output data by ordinary kriging, we derive the maps of the estimated factors and that of the permeability estimated by 2D FA-PSO method along the investigated cross-section (Fig. 11). The Fig. 9 Convergence of the 2D PSO based factor analysis procedure in Wells 1-5, Pannonian Basin, Hungary. The variation of learning parameters (c 1 plotted with dashed line, c 2 plotted with continuous line) during the iterative procedure (left panel), the decrease of data distance (energy) vs. iteration steps (right panel) Fig. 10 Regression relationship found between the first factor and the decimal logarithm of permeability measured in mD in Wells 1-5, Pannonian Basin, Hungary   Fig. 11 Result of permeability prediction using 2D PSO assisted factor analysis for multiple wells. Top panel shows the section of the extracted first factor, the middle panel includes that of the second factor, bottom panel illustrates the estimated logarithmic permeability (originally given in mD) in Wells 1-5, Pannonian Basin, Hungary depth of potential reservoirs and sealing units can be recognized well on the image (note that the depth scale represents relative depth coordinates). Of course, care must be taken, the farther the wells are from each other, the more uncertainty there is in permeability prediction because of lithological and structural variations. The reliability of the results can be ensured with the integration of core measurements and surface geophysical measurements. The CPU time of the FA-PSO procedure for automatically determined hyperparameters is~100 h, while for fixed hyperparameters is~20 min using an octa-core processor workstation.

Calibration by core data
The results of factor analysis highly depend not just on the way of optimization but also what type and quality of the input data we have in the phase of regression analysis (Fig. 3). We test the feasibility of the permeability prediction method using the combination of wireline logs and core plug data. Let us choose a well log suite from Well-3 in the investigated area, which include the natural gamma-ray intensity (GR), potassium-thorium product (K × TH), spontaneous potential (SP), deep resistivity (R d ), shallow resistivity (R s ), bulk density (ρ b ), neutron-porosity (Φ N ) and acoustic interval time ( t) logs. From core plugs, the equivalent liquid permeability (K L ) is available, which is derived from the measured gas permeability (K G ) as where κ is a constant for a particular gas in a given rock type and P m is the mean pressure (Klinkenberg 1941). We emphasize that that laboratory measurements on core plugs and wireline logging are carried out under different conditions (e.g., temperature and pressure etc.) yielding sometimes considerable differences in predicted permeability values. Moreover, core data provide only local information about permeability, while well logs are influenced by a larger extent of the rock formation, therefore the misfit between the well log derived permeability and that of core data sometimes can be significant. To improve the strength of correlation for the variables of Eq. (12), we let the most of the observed information concentrate into only one factor. We conduct the regression analysis to find the relationship between the first factor and lab permeability on those locations where core derived permeability data is available. In the investigated depth interval, core data is available at 54 non-equidistantly sampled locations, therefore wireline logging data is gathered only from these locations and standardized to serve as the input of factor analysis. The processed interval is 18 m long. First, factor loadings are estimated by Eq. (4) only for one factor. We show the calculated and rotated factor loadings in Table 3. Based on the table, we can conclude that the extracted factor from the wireline logging dataset shows a strong correlation with the GR, K × TH, ρ b and Φ N logs, separately.
Once the factor loadings are calculated, the factor scores are optimized by PSO in the same manner as detailed in . First, we solve Eq. (5) for estimating the initial values of factor scores to define the search space of the FA-PSO process [-3, 3]. We extract only one factor from the well logging dataset, therefore the number of  (8)(9) to minimize the objective function in 200 iteration steps earlier defined in Eq. (7). We let the control parameters (c 1 , c 2 ) of the PSO algorithm change automatically. At the last iteration step, the resultant values of learning parameters are c 1 = 1.98 and c 2 = 1.91, while the data distance reaches 0.66 (Fig. 12).
To interpret the regression relationship between the first factor and the decimal logarithm of absolute permeability, first we scale the former into the range of 0 to 1. As a reference value, the equivalent liquid permeability data is used. The non-linear relationship according to the model suggested in Eq. (12) between the first scaled factor and permeability is plotted in Fig. 13. The regression coefficients are found to be a = 2.95, b = 2.10, d = 0.17 while c is fixed at 1. The Spearman's rank correlation Regression relationship found between the first factor and the decimal logarithm of core permeability measured in mD in Well-3, Pannonian Basin, Hungary coefficient is -0.78, which indicates a strong non-linear relationship. Because of having some outlying values, we can see only a moderate deviation between the well log analysis derived permeability and the data available from core plugs, nevertheless, it is well described by Eq. (12). The input (observed) well logs of factor analysis are plotted in the first eight tracks in Fig. 14, the scaled first factor log is available in track 9 of the same figure. The last track contains the decimal logarithm of permeability derived by the FA-PSO method (lg(K)_FA-PSO) and the core data (lg(K)_CORE) (note that the depth scale represents relative depth coordinates). We conclude that the two independent methods determine similar permeability distribution with depth along the processed interval. Fig. 14 The result of factor analysis in Well-3, Pannonian Basin, Hungary. The input logs of FA-PSO procedure are plotted in tracks 1-8: natural gamma-ray intensity (GR), potassiumthorium product (K× TH), spontaneous potential (SP), deep resistivity (R d ), shallow resistivity (R s ), bulk density (ρ b ), neutronporosity (Φ N ) and acoustic (P-wave) transit time (Δ t ). The scaled first factor log (F 1 ) is plotted in track 9. In track 10, lg(K)_FA-PSO represents the factor analysis derived logarithmic permeability (originally given in mD), while lg(K)_CORE shows the core derived value for the same quantity discretization scheme the so-called interval inversion method with a great overdetermination (data-to-unknowns ratio) may significantly reduce the uncertainty of the inverted parameters (Dobróka et al. 2016). Global optimization tools like PSO, SA and genetic algorithms can be effectively combined with linearized inversion techniques to solve such an inverse problem. By the successive application of the PSO technique and the interval inversion method, not only the computer processing time of PSO inversion can be reduced, but there is a possibility to calculate the estimation errors of petrophysical properties from the Jacobi matrix generated in a given linear inversion step (Abordán and Szabó 2020;Szabó and Dobróka 2020).
The problem with the estimation of permeability is that the permeability is not included at all in the theoretical tool response functions connecting the model parameters with the measured data (Alberty and Hashmy 1984). Thus, the quality check of this parameter is not possible within the inversion procedure. In the practice, the permeability is derived from the inversion results (or from different other sources), by using some empirical relation. In this study, the suggested method for permeability prediction is based on factor analysis and its results are all set in a deterministic framework. Our approach provides a quite straightforward estimation for permeability using the factor scores. However, it is possible to test the reliability of the regression model by performing noise sensitivity tests using real data. Consider the well logs and core data presented in Sect. 3.4. First we add 5% Gaussian distributed noise to the input logs, then we run the FA-PSO procedure, and we repeat this process 30 times (in different runs, the noise was always regenerated and added to the observed data). The histograms of the resultant regression coefficients (a, b, d) are plotted in Fig. 15. The average and ranges of their estimated values are: a = 2.8 ± 0.7, b = 2.25 ± 0.75, d = 0.25 ± 0.75 (coefficient c is fixed as earlier). The result of regression analyses can be seen in the bottom right panel of the figure, which shows the consistency of regression model (12). The uncertainty of permeability prediction is acceptably low and is comparable with the level of data noise. This is also proven in Fig. 16, where the error intervals of the input logs, and those of the first factor and permeability are plotted, respectively. It can be seen that the uncertainty of permeability prediction is higher in formations with lower permeability, but it is favorable that the same parameter in the hydrocarbon bearing zones can be more reliably produced. The results are also confirmed by the core data available in the processed interval.  . 16 The reliability of factor analysis based permeability estimation method tested by 30 repeated runs in Well-3, Pannonian Basin, Hungary. The input logs and their assumed error ranges (at 5% Gaussian noise) are plotted in tracks 1-8: natural gamma-ray intensity (GR), potassium-thorium product (K×TH), spontaneous potential (SP), deep resistivity (R d ), shallow resistivity (R s ), bulk density (ρ b ), neutron-porosity (Φ N ) and acoustic (P-wave) transit time (Δ t ). The scaled first factor log (F 1 ) including its possible upper and lower boundaries is plotted in track 9. In track 10, lg(K)_FA-PSO represents the factor analysis derived logarithmic permeability and its uncertainty interval, while lg(K)_CORE shows the core derived value for the same parameter

Conclusions
In the paper, we introduce a particle swarm optimization assisted factor analysis with automatically optimized control parameters to permeability prediction in hydrocarbon reservoirs. The hyperparameter estimation based PSO process requires no prior knowledge on its learning parameters. A metaheuristic search is guided for estimating both the optimal set of cognitive and social constants and the factor scores. Based on the above, we perform exploratory factor analysis by a new alternative optimization tool. The advantages of the proposed method against traditional approaches are: (I) derivative-free solution is given, (II) optimal misfit is achieved between the measurements and predictions (as global optimization is made), (III) solution is practically independent from the initial selection of factor scores, (IV) highly adaptive procedure applied to different well log suites (see Table 1-3), (V) more reliable solution is given by automatically estimating the learning parameters of PSO than those treating arbitrarily chosen control parameters, (VI) robust solution can be assured using an appropriate choose of the misfit function. More robust norms like a most frequent value based or L 1 -norm can be substituted into the Euclidean norm based objective functions . At the moment, the only serious drawback of the hyperparameter estimation method is its long computational running time. For automatically calculated hyperparameters, the CPU times increase approximately 300 times compared to those of a single run using fixed values of c 1 and c 2 . However, based on preliminary program runs, parameters c 1 and c 2 can be determined for the investigated area and then fixed for factor analysis, which may reduce the running time by orders of magnitude. In the knowledge of proper initial values, instead of the SA algorithm, a linearized optimization procedure may be proposed to reduce the running time greatly.
We have established a well usable regression model between the first factor estimated by the particle swarm optimization based factor analysis and the permeability of hydrocarbon-bearing formations. It is noted that pilot tests confirmed the validity of the same regression relation in some hydrocarbon wells in the Carpathian Foredeep, South Poland (and hydraulic conductivity in Hungary and the USA). This allows the estimation of permeability from an independent well-log-analysis method, which is based on the comprehensive interpretation of multiple wireline logs. We have tested the developed method in single wells in Hungarian hydrocarbon-bearing formations, which was extended successfully to two-dimensional application to simultaneously evaluate the data measured in multiple boreholes. The factor analysis derived permeability is compared and verified both by deterministic modeling of well logs and laboratory based measurements. As a next research step, the inertia weight w might also be automatically selected for PSO, and the factor loadings will be used as unknowns as we have done it in our earlier applications. By using fast computers and the tools of parallel computing, a generalized FA-PSO method can be developed for a reliable and fast 3D permeability prediction in different oilfields.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.