A Fully-Self-Adaptive Harmony Search GMDH-Type Neural Network Algorithm to Estimate Shear-Wave Velocity in Porous Media

Abstract

Shear wave velocity (V_S) is one of the most important parameters in deep and surface studies and the estimation of geotechnical design parameters. This parameter is widely utilized to determine permeability and porosity, lithology, rock mechanical parameters, and fracture assessment. However, measuring this important parameter is either impossible or difficult due to the challenges related to horizontal and deviation wells or the difficulty in reaching cores. Artificial Intelligence (AI) techniques, especially Machine Learning (ML), have emerged as efficient approaches for dealing with such challenges. Therefore, considering the advantage of the ML, the current research proposes a novel Fully-Self-Adaptive Harmony Search—Group Method of Data Handling (GMDH)-type neural network, named FSHS-GMDH, to estimate the V_S parameter. In this way, the Harmony Memory Consideration Rate (HMCR) and Pitch Adjustment Rate (PAR) parameters are calculated automatically. A novel method is also introduced to adjust the value of the Bandwidth (BW) parameter based on the cosine wave and each decision variable values. In addition, a variable-size harmony memory is proposed to enhance both the diversification and intensification. Our proposed FSHS-GMDH algorithm quickly explores the problem space and exploits the best regions at the late iterations. This algorithm allows for the training of the prediction model based on the P-wave velocity (V_P) and the bulk density of rock (RHOB). Applying the proposed algorithm to a carbonate petroleum reservoir in the Persian Gulf demonstrates that it is capable of accurately estimating the V_S parameter better than state-of-the-art machine learning methods in terms of the coefficient of determination (R²), Mean Square Error (MSE), and Root Mean Square Error (RMSE).

1. Introduction

Seismic velocities, including P-wave velocity (V_P) and S-wave data (V_S), play a vital role in subsurface studies. For example, geomechanical modeling or rock physical studies requires the V_S parameter to be estimated as a decisive parameter. Reservoir characterization is one of the most important parts of surveys in oil industry, and research on the V_S parameter can greatly help to investigate the precise reservoir characterization. Lithofacies and fluids’ properties are two significant and useful parameters in hydrocarbon reservoir monitoring. Previous studies have shown that shear wave velocity can help determine such parameters [1,2,3,4]. Besides, the V_S’ porosity relationship can be helpful in assessing the sorting and volume of cementation, especially for sandstone reservoirs. In surface geophysical investigations, the V_S determines the size and geometry of the subsurface structures. Therefore, measuring the V_S is essential, and several methods have been introduced to estimate this parameter. Two main methods for measuring V_S are reported as direct and indirect methods, each having advantages and drawbacks. The laboratory methods directly measure the shear wave velocity by coring samples [5]. Moreover, a well-logging instrument outlines the Dipole Sonic Imager (DSI) tool, known as one of the most common ways of measuring S-wave velocity. Unfortunately, the V_S cannot be measured using coring and well-logging operations, due to a number of problems, as follows: first, in the horizontal and deviation wells and in weak formations (e.g., sandstone formation), obtaining core samples in some of the oil and gas fields is impossible. Secondly, well-logging operating is not applicable in the same situation. Thirdly, these methods are expensive and time-consuming. Lastly, the measurement of shear wave velocity may not be recorded at certain depths, or the recorded sections may include erroneous values. Consequently, this leads to different indirect methods to be presented to address these problems.

Verifying the empirical rock physical relationship is one of the indirect approaches and the number of the relevant studies has increased [6,7,8,9,10,11,12]. Oloruntobi et al. introduced a shear-wave velocity prediction method that considered the relationship between shear wave velocity and other parameters, such as P-wave velocity (V_P) and density in sedimentary rocks [13]. A P-wave is a compressional wave and is one of the two main types of elastic body waves, called seismic waves in seismology. These waves are also primary waves because they are the first waves to be recorded by seismographs and devices. This type of wave is used in subsurface studies to detect the rock, fluid, pressure, and overpressure. Maleki et al. employed correlations and artificial intelligence methods to estimate a shear wave velocity from the petrophysical logs in a well drilled in carbonate formation [14]. The relationships provided by these methods were usually specific to the area or formation, and the shear wave velocity could be estimated by localizing their coefficients in some cases.

Geostatistical analysis, to which more attention has yet to be paid, can be useful to estimate the V_S. This method establishes a spatial correlation between the variables, normalizes the database, and modifies the estimation equations [15]. Thanks to immanence advances in science and technology, Artificial Intelligent (AI) techniques have been used to estimate V_S. Recently, ML methods have become popular and been widely used, because they are inexpensive, fast track, and practical without restrictions concerning the horizontal or deviation position of the wells [16]. These methods have been used in many studies to predict the shear wave velocity by employing a different range of artificial intelligent methods, such as Artificial Neural Networks (ANNs) [17,18], Support Vector Regression (SVR) and neural networks [19,20,21], soft computing [22], machine learning, and colony-fuzzy inference systems [23,24]. The results of these studies show that using AI methods to predict shear wave velocity is reliable and favorable.

The present study introduces a novel hybrid method based on the Harmony Search algorithm [25] as a metaheuristic algorithm and GMDH-type neural networks to predict shear wave velocity using well-logging data. In contrast to earlier findings, this study calculates V_S by utilizing only two parameters, including V_P and density (RHOB). In this context, a Full-Self-Adaptive Harmony Search algorithm (FSHS) is designed and utilized to optimize a GMDH-type neural network. Several modifications are made in the FSHS algorithm. First, a two-dimensional Dynamic Harmony Memory Size (DHMS) strategy is proposed to control diversity during the search process. Therefore, the FSHS algorithm quickly explores the problem space and exploits the best regions at the late iterations. Secondly, all of the parameters values are adaptive and calculated automatically in the FSHS algorithm. Thirdly, a novel method is introduced to adjust the value of the BW parameter based on the cosine wave and each decision variable value.

The rest of this paper is organized as follows: The experimental methods are explained in Section 2. The materials and preliminaries of the proposed algorithm will be described in Section 3. In Section 4, the proposed hybrid algorithm is presented. The experimental results and conclusion will be discussed in Section 5 and Section 6, respectively.

2. Experimental Methods

There is a varied range of proposed experimental methods to estimate shear-wave velocity. This study uses some of methods to calculate shear wave velocity, V_S, by employing V_P and lithology. For instance, Castagna et al. [8] introduced an equation to estimate V_S as follows:

$${\mathrm{V}}_{\mathrm{S}}=0.86{\mathrm{V}}_{\mathrm{P}}-1170$$

(1)

In addition, Greenberg and Castagna introduced one of the most popular methods to estimate V_S, by considering V_P and lithology [6] (Equation (2)):

$${\mathrm{V}}_{\mathrm{S}}=\frac{1}{2}\left\{\left[{\displaystyle \sum _{\mathrm{i}=1}^{\ell }{\mathrm{X}}_{\mathrm{i}}}\left({\displaystyle \sum _{\mathrm{j}=0}^{{\mathrm{N}}_{\mathrm{i}}}{\mathrm{a}}_{\mathrm{i}\mathrm{j}}{\mathrm{V}}_{\mathrm{p}}^{\mathrm{j}}}\right)\right]+{\left[{\displaystyle \sum _{\mathrm{i}=1}^{\ell }{\mathrm{X}}_{\mathrm{i}}}{\left({\displaystyle \sum _{\mathrm{j}=0}^{{\mathrm{N}}_{\mathrm{i}}}{\mathrm{a}}_{\mathrm{i}\mathrm{j}}{\mathrm{V}}_{\mathrm{p}}^{\mathrm{j}}}\right)}^{-1}\right]}^{-1}\right\}$$

(2)

where $\ell$ refers to the number of pure components in terms of lithology; X_i presents volume proportion of lithological constituents; a_ij is empirical regression coefficients relying upon lithology (

Table 1. Regression coefficients for V_S prediction for the Greenberg—Castagna relations.

Lithology	ai2	ai1	ai0
Sandstone	0	0.80416	−0.85588
Limestone	−0.05508	1.01677	−1.03049
Dolomite	0	0.58321	−0.07775
Shale	0	0.76969	−0.86735

); N_i refers to the order of polynomial for constituent I; ${\mathrm{V}}_{\mathrm{p}}^{\mathrm{j}}$ is the water-saturated P-wave velocity in the j rock facies. Finally, V_S is S–wave velocities (km/s) in composite multi-mineral brine-saturated rock.

Castagna et al. [8] provided three other equations based on V_P and lithology (Equations (3)–(5)):

$${\mathrm{V}}_{\mathrm{S}}=0.5508{\mathrm{V}}_{\mathrm{P}}^2+1.0116{\mathrm{V}}_{\mathrm{P}}-1.0305$$

(3)

$${\mathrm{V}}_{\mathrm{S}}=0.5832{\mathrm{V}}_{\mathrm{P}}-0.0777$$

(4)

$${\mathrm{V}}_{\mathrm{S}}={0.77\mathrm{V}}_{\mathrm{P}}-0.8674$$

(5)

These equations are used for a specific purpose. For instance, Equation (3) is used for limestone, and Equations (4) and (5) are used for dolomite and shale, respectively. Furthermore, Brocher et al. [7] introduced a model based on V_P to calculate the shear wave velocity for each formation (Equation (6)):

$${\mathrm{V}}_{\mathrm{S}}=0.7858-1.2344{\mathrm{V}}_{\mathrm{P}}^2-0.1238{\mathrm{V}}_{\mathrm{P}}^3+0.0006{\mathrm{V}}_{\mathrm{P}}^4$$

(6)

3. Materials

3.1. Harmony Search (HS) Algorithm

Traditional optimization algorithms, such as linear programming (LP), non-linear programming (NLP), and dynamic programming (DP), have taken on major roles in solving optimization problems. Nevertheless, their drawbacks generate a demand for other types of algorithms, such as meta-heuristic/heuristic optimization approaches. These algorithms have been identified as efficient, intelligent techniques to solve optimization issues. In recent years, various metaheuristic algorithms (MAs) have been proposed to solve a wide range of complex real-world problems [26,27,28]. The harmony search (HS) algorithm is a type of metaheuristic optimization algorithm, that was designed by Geem et al. [25]. Inspired by the improvisation of music players, they first play music randomly with existing instruments. This harmony is kept in the musician’s memory. In the next stage, according to the harmony in their memory, the musician plays new music that has changed from the previous one. Generally, the HS algorithm consists of three steps: (a) harmony memory initialization; (b) improvisation; and (c) updating the harmony memory. The harmony memory initialization step produces the initial harmonies (solutions), randomly. Then, in step (b), an iterative process is completed to enhance the quality of solutions, by making new solutions according to three operators, including harmony memory consideration, pitch adjustment, and random consideration. Finally, in step (c), the best-found harmony is selected based on objective function(s) as the optimum solution.

The HS algorithm has been applied to solve a wide range of real world optimization problems. In [29], a modified harmony search algorithm was introduced and applied to overcome the problem of combined heat and power economic dispatch. In order to solve the Team Orienteering Problem (TOP), a new version of HS algorithm was proposed in [30]. In this way, the authors proposed an innovative method called the Similarity Hybrid Harmony Search (SHHS) algorithm. A variant of the original harmony search algorithm was developed in [31] for the optimum design of steel frames. The authors applied the HS algorithm to find minimum weight frames by choosing appropriate sections from a standard steel sections set. In [32], a discrete version of HS algorithm, named DHS, was introduced to solve a Flexible Job-shop Scheduling Problem (FJSP). The proposed DHS algorithm was utilized to minimize a weighted combination of the maximum completion time and the earliness mean. A hybrid variant of the Harmony Search (HHS) algorithm was introduced by the authors of [33]. The authors hybridized the neighborhood search and the global-best search mechanism of the Particle Swarm Optimization (PSO) algorithm with the HS algorithm to solve the problem of discrete truss-structure optimization. An adaptive harmony search algorithm was designed to achieve an optimal virtual machine (VM) consolidation model that decreased the power consumption of the datacenter in [34]. In order to obtain the optimum solution for the proportional-integral-derivative (PID) controller, M’zough et al. proposed a Self-Adaptive Global Harmony Search (SGHS) algorithm [35]. Sasha et al. [36] introduced a modified HS algorithm called Supervised Filter Harmony Search Algorithm (SFHSA), to eliminate redundant/irrelevant features without any significant effect on the process of classification for a Facial Emotion Recognition (FER) task.

3.2. Group Method of Data Handling (GMDH)

Introduced by Alexey G. Ivakhnenko [37], the GMDH-type neural network is a self-organized method that combines both self-organizing and multilayer neural network benefits to fit the data [38]. The GMDH tries to automatically identify the best model for structuring the network and to integrate inductive techniques, biological neuron techniques, probability, and black-box concepts [39]. Therefore, according to the GMDH, a model is a set of neurons made by pairs based on the quadratic polynomial combination of inputs or previous layers’ outputs. The main goal of this method is to find function f to estimate in a feed-forward network according to a set of N inputs X = {x₁, x₂,…, x_N} and one output Y, based on a quadratic polynomial transfer function. The Volterra series is a non-linear transfer function used in GMDH neural networks, shown in Equation (7). In this equation, the regression method is used to calculate the coefficients a_i for pairs of input variables (x_i, x_j) [39]. As an example, a two-variable second-degree polynomial, derived from the Volterra series, is presented in Equation (8):

$$\widehat{\mathrm{Y}}={\mathrm{a}}_0+{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{a}}_{\mathrm{i}}{\mathrm{x}}_{\mathrm{i}}+}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{N}}{\mathrm{a}}_{\mathrm{i}\mathrm{j}}{\mathrm{x}}_{\mathrm{i}}{\mathrm{x}}_{\mathrm{j}}+}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{N}}{\displaystyle \sum _{\mathrm{k}=1}^{\mathrm{N}}{\mathrm{a}}_{\mathrm{i}\mathrm{j}\mathrm{k}}{\mathrm{x}}_{\mathrm{i}}{\mathrm{x}}_{\mathrm{j}}{\mathrm{x}}_{\mathrm{k}}+}}}\dots$$

(7)

$$\mathrm{f}({\mathrm{x}}_{\mathrm{i}},{\mathrm{x}}_{\mathrm{j}})={\mathrm{a}}_0+{\mathrm{a}}_1{\mathrm{x}}_{\mathrm{i}}+{\mathrm{a}}_2{\mathrm{x}}_{\mathrm{j}}+{\mathrm{a}}_3{\mathrm{x}}_{\mathrm{i}}^2+{\mathrm{a}}_4{\mathrm{x}}_{\mathrm{j}}^2+{\mathrm{a}}_5{\mathrm{x}}_{\mathrm{i}}{\mathrm{x}}_{\mathrm{j}}$$

(8)

4. Methods

In this section, the proposed hybrid FSHS-GMDH algorithm is explained in detail. First, the Fully-Self-Adaptive Harmony Search (FSHS) algorithm is described step by step. Then, we will explain the hybridization scheme.

4.1. Fully-Self-Adaptive Harmony Search (FSHS) Algorithm

In the FSHS algorithm, several modifications have been made to the original HS algorithm. In this way, a new method is presented, based on the Harmony Memory (HM) with variable size and other values of decision variables. The FSHS algorithm aims to explore the problem space quickly and exploits the best regions at the late iterations. The values of the Harmony Memory Consideration Rate (HMCR) and the Pitch Adjustment Rate (PAR) parameters are also automatically calculated. In addition, a novel method is introduced to adjust the value of the Bandwidth (BW) parameter based on the cosine wave [26] and each decision variable value. The flowchart of the novel Full-Self-Adaptive Harmony Search algorithm is illustrated in Figure 1.

4.1.1. Variable-Size Harmony Memory

In the main HS, the harmony memory size is considered constant. It leads to the use of all of the memory solutions in producing new harmonies. Although this increases the diversity in the production of new harmonies in the initial iterations, it slows down the convergence rate in the final iterations [40].

We used a variable-size memory in the proposed method to address the problem, as shown in Equation (9):

$$\mathrm{H}\mathrm{M}{\mathrm{S}}_{\mathrm{g}}=\lceil \left(1-{\left(\frac{0.9}{\mathrm{N}\mathrm{I}}\times \mathrm{g}\right)}^2\right)\times \mathrm{H}\mathrm{M}{\mathrm{S}}_{\max }\rceil$$

(9)

where HMS_g is the harmony memory size in generation g; HMS_max is the maximum size of harmony memory; NI is the number of iterations in the algorithm; and g is the desired generation.

As all of the memory harmonies are used at the beginning of the search and during the searches’ progress, the number of harmonies will decrease gradually. Lower-quality harmonies are not physically removed from the harmony memory and are utilized in the random selection step. In the initial iterations of this strategy, using all of the harmonies in the memory leads to maintenance of the ability of the algorithm to keep the diversity by generating various harmonies. Besides, the intensification increases in the final iterations by focusing on the best harmonies. Moreover, as a result of increasing the diversity in the production of a new solution and increasing the algorithm’s capability in the global search in the proposed method, instead of using the values of the same decision variable, all of the decision variables values are used after selecting a harmony from the memory. A similar method is used in GHS [41] and IGHS [42]. However, this process is completed in the Pitch Adjustment section and only the values of the decision variables of the best harmony are used. Although this method greatly increases the convergence rate in the initial iterations, it reduces the exploitation ability of the algorithm from the middle iterations onwards.

Hence, if the optimal values of the variables are not a close decision, selecting the values of other variables prevents convergence in the final iterations [42]. Therefore, a reduction mechanism is utilized in the FSHS to use the values of the other decision variables for solving the problems (Equations (10)–(12)):

$$\mathsf{\mu }=\lceil \left(1-\left(\frac{0.9}{\mathrm{N}\mathrm{I}}\times \mathrm{g}\right)\right)\times \mathrm{N}\rceil$$

(10)

$${\mathrm{S}}^{\mathsf{\mu }}\subset {\mathrm{X}}^{\mathrm{i}},\mathrm{where}{\mathrm{X}}^{\mathrm{i}}=\left\{{\mathrm{x}}_1^{\mathrm{i}},{\mathrm{x}}_2^{\mathrm{i}},{\mathrm{x}}_3^{\mathrm{i}},\dots ,{\mathrm{x}}_{\mathrm{N}}^{\mathrm{i}}\right\}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{i}=\mathrm{U}(1,\mathrm{H}\mathrm{M}{\mathrm{S}}_{\mathrm{g}})$$

(11)

$${\mathrm{X}}_{\mathrm{j}}^{\mathrm{n}\mathrm{e}\mathrm{w}}={\mathrm{X}}_{\mathrm{k}}^{\mathrm{i}},\mathrm{where}{\mathrm{X}}_{\mathrm{k}}^{\mathrm{i}}\in {\mathrm{S}}^{\mathsf{\mu }}\mathrm{and}\mathrm{k}=\mathrm{U}(1,\mathsf{\mu })$$

(12)

where ${\mathrm{S}}^{\mathsf{\mu }}$ is a randomly selected of the harmony decision variables set which includes μ members; Xⁱ and k is an integer random number between 1 and μ; and the decision variable ${\mathrm{X}}_{\mathrm{k}}^{\mathrm{i}}$ is a member of the ${\mathrm{S}}^{\mathsf{\mu }}$. In other words, there is a variety of choice for the values in the initial iterations. From the middle iterations onwards, the focus is on the values of the same decision variables.

4.1.2. Harmony Memory Consideration Rate (HMCR) and Pitch Adjustment Rate (PAR)

The two parameters, HMCR and PAR, have an effective role in the efficiency of the HS algorithm and extensive research has been completed on the automatic adjustment of these parameters. Mahdavi et al. [43] introduced a new mechanism for setting the PAR parameter. In their method, the PAR values are calculated between PAR_min and PAR_max for each generation of harmonies. The drawback of this method is that two new parameters need to be set. The proposed FSHS approach uses the same method (Equation (13)), except that the PAR_min and PAR_max values are omitted. Besides, Equation (14) is used to set the HMCR parameter:

$$\mathrm{P}\mathrm{A}{\mathrm{R}}_{\mathrm{g}}=\left(\frac{1}{\mathrm{N}\mathrm{I}}\times \mathrm{g}\right)$$

(13)

$$\mathrm{H}\mathrm{M}\mathrm{C}{\mathrm{R}}_{\mathrm{g}}=\left(\frac{1}{\mathrm{N}\mathrm{I}}\times \mathrm{g}\right)$$

(14)

4.1.3. Bandwidth (BW)

The studies show that determining the appropriate value for the bandwidth, BW, greatly affects the performance of the HS algorithm [19]. To be more accurate, a large BW in the initial iterations can significantly help to search the problem space and find better areas. In addition, small BW in the final iterations contributes to a local search in the neighborhood of the solutions found by the searching algorithm. So far, many methods have been proposed to adjust this parameter; however, the proper adjustment of this parameter to solve various problems is much more difficult than other parameters. In this research, a new mechanism for adjusting the BW parameter is presented, based on the proposed method in the SCA [26] algorithm to overcome such challenges. Instead of using a constant value for the BW, a variable domain is used in the proposed approach. The idea behind this method is to increase the probability of using better values by making alternating changes in the base BW range. For this purpose, the difference between the minimum and maximum values of the memory harmony decision variables is used as the value of the base BW. Then, the BW value is changed periodically, using the trigonometric function cosine and Equations (15)–(18). Furthermore, each decision variable has a specific BW in this method:

$$\overrightarrow{\mathsf{\beta }}=\left(-\left(\frac{\mathrm{N}\mathrm{I}}{10}\mathsf{\pi }\right),\left(\frac{\mathrm{N}\mathrm{I}}{10}\mathsf{\pi }\right)\right)$$

(15)

$${\mathrm{D}}_{\mathrm{j}}=\mathrm{M}\mathrm{a}\mathrm{x}({\mathrm{X}}_{\mathrm{j}})-\mathrm{M}\mathrm{i}\mathrm{n}({\mathrm{X}}_{\mathrm{j}})$$

(16)

$$\mathrm{B}{\mathrm{W}}_{\mathrm{j}}^{\mathrm{g}}={\mathrm{D}}_{\mathrm{j}}+({\mathrm{D}}_{\mathrm{j}}\times \cos ({\mathsf{\beta }}_{\mathrm{g}}))$$

(17)

$${\mathrm{X}}_{\mathrm{j}}^{\mathrm{n}\mathrm{e}\mathrm{w}}={\mathrm{X}}_{\mathrm{j}}^{\mathrm{n}\mathrm{e}\mathrm{w}}+\mathrm{U}(-1,+1)\times \mathrm{B}{\mathrm{W}}_{\mathrm{j}}^{\mathrm{g}}$$

(18)

where β is a vector with length NI between $-(\frac{\mathrm{NI}}{10}\mathsf{\pi })$ and $+(\frac{\mathrm{NI}}{10}\mathsf{\pi })$ with intervals $(\frac{2\times \mathrm{NI}}{10\times \mathrm{NI}}\mathsf{\pi })$; D_j is the difference between the maximum and minimum values of X_j in the harmony memory; and $\mathrm{B}{\mathrm{W}}_{\mathrm{j}}^{\mathrm{g}}$ is the bandwidth of the decision variable X_j in the generation of g. In addition, ${\mathrm{X}}_{\mathrm{j}}^{\mathrm{n}\mathrm{e}\mathrm{w}}$ is the jth decision variable of the new harmony ${\mathrm{X}}^{\mathrm{n}\mathrm{e}\mathrm{w}}$ and U(−1, +1) is a random value with uniform distribution between −1 and +1.

4.1.4. Random Selection

The proposed FSHS algorithm, adopted from [40], utilizes a similar approach for the random selection step. In doing so, instead of generating random values between the upper bound (UB) and lower bound (LB) in this step, we utilize the mean value of the decision variable and a random value between the maximum and minimum values of the decision variables of the HM (Equation (19))

$${\mathrm{X}}_{\mathrm{j}}^{\mathrm{n}\mathrm{e}\mathrm{w}}=\overline{{\mathrm{X}}_{\mathrm{j}}}\pm \mathrm{U}(-1,+1)\times \left(\frac{\mathrm{M}\mathrm{a}\mathrm{x}({\mathrm{X}}_{\mathrm{j}})-\mathrm{M}\mathrm{i}\mathrm{n}({\mathrm{X}}_{\mathrm{j}})}{2}\right)$$

(19)

where ${\mathrm{X}}_{\mathrm{j}}^{\mathrm{n}\mathrm{e}\mathrm{w}}$ is the jth decision variable of new harmony ${\mathrm{X}}^{\mathrm{n}\mathrm{e}\mathrm{w}}$ and $\overline{{\mathrm{X}}_{\mathrm{j}}}$ is the mean value of jth decision variable in HM.

4.2. Hybridization of the FSHS and GMDH-Type Neural Network Algorithms

The process of hybridization of our proposed FSHS-GMDH algorithm is explained in this section. As mentioned earlier, the GMDH-type neural network conducts a self-organized evolutionary mechanism to find an optimum structure, layer by layer. Therefore, in the proposed hybrid FSHS-GMDH algorithm, the FSHS algorithm would be launched to optimize each layer, and the coefficients of all the neurons in a quadratic polynomial equation form are considered as decision variables and fine-tuned by the FSHS algorithm. In addition, we consider the Mean Square Error (MSE) as the cost function that must be minimized. A schematic view of the proposed FSHS-GMDH algorithm is shown in Figure 2.

5. Results and Discussion

In this section, the implementation condition, dataset, and inputs are explained. After that, the results will be discussed in detail. All of the algorithms and methods are implemented using MATLAB software (version R2017b). Furthermore, the following configurations and settings are considered in our study of the algorithms. For the FSHS-GMDH algorithm and the original GMDH, each maximum number of the layers and neurons is set to 4 and 10, respectively. In addition, the number of neurons in the hidden layer is set at 8 (4-8-1) for the multilayer perceptron (MLP) model. Furthermore, R², MSE, and Root Mean Squared Error (RMSE) are calculated, using Equations (20)–(22) to evaluate and compare the proposed hybrid algorithm and other methods. We used only two parameters, including density (RHOB) and V_P in this study to estimate the shear wave velocity. Figure 3 shows the inputs data and original shear wave velocity. Equations (20)–(22) are as follows:

$${\mathrm{R}}^2=\sqrt{\frac{\mathrm{N}\times {\displaystyle \sum _{\mathrm{t}=1}^{\mathrm{N}}({\mathrm{X}}_{\mathrm{t}}\times {\mathrm{Y}}_{\mathrm{t}})}-{\displaystyle \sum _{\mathrm{t}=1}^{\mathrm{N}}({\mathrm{X}}_{\mathrm{t}})}\times {\displaystyle \sum _{\mathrm{t}=1}^{\mathrm{N}}({\mathrm{Y}}_{\mathrm{t}})}}{\sqrt{[\mathrm{N}\times {\displaystyle \sum _{\mathrm{t}=1}^{\mathrm{N}}({\mathrm{X}}_{\mathrm{t}}^2)}-({\displaystyle \sum _{\mathrm{t}=1}^{\mathrm{N}}{({\mathrm{X}}_{\mathrm{t}})}^2}]\times [\mathrm{N}\times {\displaystyle \sum _{\mathrm{t}=1}^{\mathrm{N}}({\mathrm{Y}}_{\mathrm{t}}^2)}-({\displaystyle \sum _{\mathrm{t}=1}^{\mathrm{N}}{({\mathrm{Y}}_{\mathrm{t}})}^2}]}}}$$

(20)

$$\mathrm{M}\mathrm{S}\mathrm{E}=\frac{{\displaystyle \sum _{\mathrm{t}=1}^{\mathrm{N}}{({\tilde{\mathrm{Y}}}_{\mathrm{t}}-{\mathrm{Y}}_{\mathrm{t}})}^2}}{\mathrm{N}}$$

(21)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{\mathrm{N}}\times {\displaystyle \sum _{\mathrm{t}=1}^{\mathrm{N}}{({\tilde{\mathrm{Y}}}_{\mathrm{t}}-{\mathrm{Y}}_{\mathrm{t}})}^2}}$$

(22)

In our research, the used shear wave velocity is original and measured by the dipole sonic log in the well-logging operation. The correlation between the mentioned parameters and the shear wave velocity must be evaluated in the first step. Figure 4 illustrates that the correlation between the V_S and two selected input data are 0.71 and 0.86 for RHOB and V_P, respectively. Therefore, they can be reliable as the input data.

5.1. Estimation Shear Wave Velocity Comparison

Figure 4 illustrates the plot diagrams for the FSHS-GMDH, MLP, and original GMDH models. The vertical axis represents the estimated values as the output value of the model and the horizontal axis shows the target values as the actual observations. Figure 5a–c illustrate the results of the proposed FSHS-GMDH. Figure 5d–f represents the MLP model results, and the original GMDH model results are shown in Figure 5g–i, respectively. According to the results, it can be seen that our proposed FSHS-GMDH algorithm outperforms the MLP and the original GMDH models in shear wave velocity estimation. In addition, the predicted results of the shear wave velocity and the error distribution for test and total data (train and test) are shown in Figure 6 and Figure 7, respectively. As can be seen, the proposed FSHS-GMDH algorithm is capable enough to estimate the shear wave velocity, according to the test data.

5.2. Comparison with Experimental Methods

The statistical analysis of the shear wave velocity is calculated by experimental methods (Equations (20)–(22)) in terms of the coefficient of determination R², which is presented in

Table 2. Statistical analysis (R²) of shear wave velocity estimation calculated by experimental methods Equations (1)–(6).

Castagna et al. [8] Equation (1)	Greenberg and Castagna [6] Equation (2)	Castagna et al. [8] Equations (3)–(5)			Brocher [7] Equation (6)	FSHS-GMDH
		Limestone	Dolomite	Shale
0.8130	0.8103	0.7962	0.8054	0.8051	0.4945	0.9688

. According to the results in this Table, it is clear that the proposed FSHS-GMDH algorithm achieves the best performance for all of the data with R² = 0.9688. In addition, the second-best result is achieved with the experimental method introduced by Castagna et al. [6]. However, there is no significant difference between the experimental methods’ results as presented in Table 2.

5.3. Comparison with Machine Learning Methods

In order to have a fair comparison, the performance of the FSHS-GMDH algorithm is investigated using the train and test data. In doing so, the MLP, GMDH, and FSHS-GMDH algorithms are firstly trained to produce the optimum models. Then, obtained from the training phase, the best setup of each algorithm is employed to estimate the shear wave velocity, according to the test data samples.

Table 3. Statistical analysis of shear wave velocity estimation calculated by ML methods.

	Train Data			Test Data			All Data
Algorithm	R2	RMSE	MSE	R2	RMSE	MSE	R2	RMSE	MSE
MLP	0.8691	93.48	8738.92	0.8512	103.82	10778.8	0.8634	96.70	9351.77
GMDH	0.8935	83.64	6995.87	0.8647	101.49	10301.39	0.8834	89.38	7988.96
FSHS-GMDH	0.9728	43.36	1880.23	0.9514	55.82	3116.92	0.9688	46.12	2127.93

represents the statistical results of shear wave velocity estimation. The results are presented in R², RMSE, and MSE for the train and test data. It can be seen from Table 3 that the proposed FSHS-GMDH algorithm provides a better performance than the MLP and the original GMDH methods. The FSHS-GMDH algorithm provides the lowest values of RMSE (43.36), MSE (1880.23), and the highest value of R² (0.9728) for the train data. Furthermore, according to the results for the test data, the best performance is achieved by the FSHS-GMDH model with the lowest values of MSE (3116.92), RMSE (55.82), and the highest value of R² (0.9514). However, the original GMDH model shows better performance for the train data; there is no significant difference between the performance of the original GMDH and MLP models for the test data.

6. Conclusions

This study estimated shear wave velocity as one of the most important parameters in reservoir geophysical characterization via the Full-Self-Adaptive Harmony Search-GMDH neural network algorithm in a carbonate petroleum reservoir in the Persian Gulf. In doing so, a novel version of the Harmony Search (HS) algorithm was introduced and applied to optimize a GMDH-type neural network. In addition, only two of the parameters, including density (RHOB) and P-wave velocity (V_P), were utilized as inputs, as they have the highest correlation with the S-wave velocity, 0.71, and 0.86, respectively. The results proved that estimating S-wave data through our novel algorithm with R² ≈ 0.97 is more accurate and cost-effective than similar studies. Previous studies used more than two inputs of data that led to spending more time and higher costs; however, their results were not as accurate as the results in this study. Nevertheless, from an optimization perspective, metaheuristic algorithms concentrate on providing a comprehensive search to find a global optimum and it may lead to the model become over-fit. In order to overcome this problem, it is crucial to conduct a deep study to enhance the generalizations in the future work.