# A Modified Whale Optimization Algorithm and Its Application in Seismic Inversion Problem

Authors: Xiaodan Liang; Siwen Xu; Yang Liu; Liling Sun
Journal: Mobile Information Systems
 (2022)

---

## Abstract

The whale optimization algorithm (WOA) is a popular swarm intelligence algorithm which simulates the hunting behavior of humpback whales. WOA has the deficiency of easily falling into the local optimal solutions. In order to overcome the weakness of the WOA, a modified variant of WOA called OCDWOA is proposed. There are four main operators introduced into the OCDWOA to enhance the search performance of WOA. The operators include opposition-based learning method, nonlinear parameter design, density peak clustering strategy, and differential evolution. The proposed algorithm is tested on 19 optimization benchmark functions and a seismic inversion problem. OCDWOA is compared with the classical WOA and three typical variants of WOA. The results demonstrate that OCDWOA outperforms the compared algorithms in terms of obtaining the global optimal solution.


---

## Body

## 1. Introduction

In recent years, the nature-inspired algorithms have attracted lots of attention from researchers. Swarm intelligence (SI) algorithm is a main kind of nature-inspired algorithms. Swarm intelligence algorithms are designed by simulating the behaviors of animals or plants. The widely used SI algorithms include particle swarm optimization (PSO) [1], Ant colony optimization (ACO) [2], artificial bee colony algorithm (ABC) [3], grey wolf optimization (GWO) [4], wolf pack algorithm (WPA) [5], brain storm optimization [6], and whale optimization algorithm (WOA) [7].The whale optimization algorithm (WOA) is proposed in 2016. WOA searches the global optimal solution by simulating the foraging behaviors of the whales, including encircling prey, researching prey, and attacking prey. Many studies about WOA have been presented. Generally, these studies can be classified into two categories. One is to apply applying the WOA to solve some practical issues. The other one is to improve the performance of WOA with some strategies. A modified WOA named chaotic WOA (CWOA) [8] was proposed by Kaur and Arora in 2017. The CWOA introduces chaos theory into WOA. The main parameter of WOA is tuned by chaotic map to control exploration and exploitation. In [9], Lévy flight trajectory-based WOA (LWOA) is presented. The Lévy flight trajectory is used to increase the diversity of the population and enhance the performance of jumping out of local optimal solutions. The balanced variant of WOA (BWOA) proposed in [10] introduced Lévy flight and chaotic local search into WOA to enhance the inclusive exploratory and neighborhood-informed capacities of the conventional technique. Chen [11] proposed a modified WOA based on adaptive convergence and Lévy features (IWOA). The adaptive convergence operator is used to balance the local and global optimization ability and the Lévy flight mechanism is introduced to improve global searching ability. For solving the large-scale global optimization problems, a modified self-adaptive WOA (MWOA) [12] is proposed. A nonlinear dynamic strategy is applied to update the control parameter and balance the exploration and exploitation abilities. The Lévy-flight strategy is utilized to enhance the search performance and a quadratic interpolation method is adopted to improve the exploitation ability of the algorithm. For solving high-dimensional function optimization problems, A WOA based on Lamarckian learning (WOALam) [13] brought in the evolutionary theory of Lamarck to enhance the local search. Trivedi et al. [14] proposed an adaptive whale optimization algorithm (AWOA) with a randomization technique. The results show that the performance of AWOA is better than the performance of the standard WOA. Sun et al. proposed a whale optimization algorithm based on quadratic interpolation (QIWOA) [15]. A modified exploration process and quadratic interpolation are used to improve the exploitation ability and the search accuracy. A reinforced variant of WOA named RDWOA [16] was proposed in 2019. A strategy of random spare or random replacement is used to enhance the convergence speed. In addition, a strategy of double adaptive weight is introduced to improve the exploration ability in the early stage and exploitative ability in the later stage. Trivedi et al. [17] proposed a hybrid PSO-WOA. The PSO operator is used in exploitation phase and the basic WOA operator is adopted to explore in uncertain environments. The experimental results reveal the effectiveness of PSO-WOA compared to standard PSO and WOA. Wei et al. [18] proposed a modified WOA based on different searching paths and perceptual disturbance (PDWOA). Several spiral curves have been tested and the results show that the equal-pitch of Archimedean spiral curve is more capable than other types of spiral curve. Mafarja and Mirjalili [19] proposed a hybrid SA-WOA algorithm. In order to enhance the exploitation of WOA, the simulated annealing (SA) algorithm is introduced into the algorithm by searching the most promising regions located by basic WOA. An improved whale optimization algorithm (IWOA) [20] is proposed in 2016. Inertia weight, as a new parameter, is introduced to adjust the best solution in each iteration. Abdel-Basset et al. [21] proposed a memetic algorithm combining the WOA with a Tabu Search for solving quadratic assignment problems.In addition, WOA has been applied to many engineering problems. Oliva et al. [22] used an improved chaotic whale optimization algorithm to estimate the parameters of solar cells. The design problem involves the solution of the complex nonlinear and multimodal objective functions. The proposed algorithm solves the design problem with the chaotic map and automatically adapt parameters. Zhang et al. [23] proposed an enhanced whale optimization algorithm (EWOA) to optimize the flow shop construction scheme. The enhanced WOA with Lévy flight is used to improve the robustness and efficiency of the existing construction stage and zone optimization approaches. Bentouati et al. [24] proposed a new variant WOA with pattern search algorithm (PS). The results prove the effectiveness and the performance of the proposed algorithm for solving the optimal power flow problem. Abd Elaziz and Oliva [25] proposed an enhanced WOA based on opposition (OBWOA). The opposition-based learning is used to enhance the exploration of the search space. OBWOA has the exploration abilities for different benchmark optimization functions and has been applied to estimate three different diode models of solar cells. Aljarah et al. [26] applied WOA to neural networks training. The proposed algorithm has been tested on 20 datasets. The results showed that the proposed method can obtain the best solutions and outperform the other compared algorithms on convergence speed. Medani et al. [27] proposed a modified WOA to solve the optimal reactive power dispatch problem. The results show that the proposed algorithm is more desirable than two comparative algorithms.Prestack seismic inversion technology is an important seismic inversion technology to solve the problems of reservoir lithology and fluid prediction in oil and gas exploration and development. Physical properties obtained by general prestack inversion methods provide many detailed subsurface characteristics, but the method is not mature and limited by the high resolution, low speed, and poor stability of prestack inversion. Currently, this technology is in the research stage. Some inevitable technology problems should be solved before its large-scale applications [28, 29].Many researches on solving the problems of prestack inversion by optimization algorithms have been studied. Xia et al. [30] and Sen [31] improved the speed and stability of prestack seismic waveform inversion by using the step inversion method. Simulated annealing (SA) and genetic algorithm (GA) have been successfully used in prestack seismic waveform inversion. Ingber and Rosen [32] conducted a comparative study on SA and GA. The experiment results proved that both algorithms could achieve global optimization and the speed of fast simulated annealing (VFSA) was the fastest. Liu et al. [33] proposed a new adaptive genetic algorithm for seismic inversion to overcome slow convergence rate and falling into local optimal solutions. Wu et al. [34] proposed a parameter inversion method based on improved differential evolution algorithm. The proposed method was effective in solving parameter inversion problems of prestack seismic data.WOA has a good convergence rate, but the performance of WOA in finding the global optimal solution of the multimodal problems with many local optimal solutions is not ideal. In order to enhance the global search ability and solve the prestack seismic waveform inversion, a modified whale optimization algorithm named OCDWOA is proposed in this paper. There are four operators used in the proposed algorithm. Nonlinear parameter design is introduced to balance the ratio of exploitation phase and exploration phase. Density peak clustering strategy is adopted to decompose population into numbers of groups for enhancing the diversity of population. Opposition-based learning method and differential evolution are used to improve the search performance. Experiments on 19 optimization benchmark functions and a seismic inversion problem indicate the effectiveness of the proposed algorithm. The proposed algorithm OCDWOA represents outstandingly the experiments. The results demonstrate that the proposed algorithm OCDWOA can solve the seismic inversion problem successfully.The main contributions are summarized as follows:(1)
The opposition-based learning method adopted by OCDWOA could avoid the solutions trapping into local optimal, which improves the ability of global search.(2)
The nonlinear parameter design ofa in OCDWOA could achieve the global exploration in early stage and local exploitation in late stage.(3)
Density peak clustering strategy is adopted to decompose population into numbers of groups, which enhances the diversity of population.(4)
The proposed OCDWOA has an outstanding performance of solving seismic inversion problem, which verifies its universality for solving practical engineering problems.The rest of this paper is organized as follows. Section2 describes the WOA algorithm. Section 3 explains the proposed algorithm and the operators. The experimental results and comparisons are shown in Section 4. Finally, in Section 5, the conclusions are presented.


## 2. Whale Optimization Algorithm

WOA is a kind of population-based meta-heuristic algorithm. WOA simulates the hunting mode of humpback whales, which includes encircling prey, bubble-net attacking method, and search for prey.


### 2.1. Encircling Prey

WOA supposes that the optimal solution with the best fitness in current population is the prey. All the other solutions will encircle the prey and update themselves toward the optimal solution. The mathematical model of encircling behavior is as follows:(1)D⇀=C⇀⋅X⇀∗t−X⇀t,(2)X⇀t+1=X⇀∗t−A⇀⋅D⇀,where t indicates the number of the current iteration, X is the current solution, X∗ is the best solution in the population, and A⇀ and C⇀ are the coefficient vectors calculated as follows:(3)A⇀=2a⇀⋅r⇀−a⇀,C⇀=2r⇀,where a⇀ is linearly decreased from 2 to 0 over the iterations and r⇀ is a random vector in [0, 1].


### 2.2. Bubble-Net Attacking Method

In order to simulate the behavior of bubble-net attacking, shrinking encircling mechanism and spiral updating position are designed as follows.


#### 2.2.1. Shrinking Encircling Mechanism

The mathematical model of shrinking encircling mechanism is defined as equations (1) and (2). Because of a⇀ is linearly decreased from 2 to 0, the value of A⇀ is in the range of [−1, 1]. The current solution can update itself toward the best solution with a random step.


#### 2.2.2. Spiral Updating Position

The new position of current solution updated with the best solution is obtained as follows:(4)D⇀′=X⇀∗t−X⇀t,(5)X⇀t+1=D⇀′⋅ebl⋅cos2πl+X∗t,where b is constant for defining logarithmic spiral and l is randomly generated number in the range [−1, 1]. To balance the two models, a selection parameter p is set. The selection parameter p in WOA is 0.5; there is a probability of 50% to choose either the shrinking encircling mechanism or the spiral updating position model. The finial mathematical model of bubble-net attacking is as follows:(6)X⇀t+1=X⇀∗t−A⇀⋅D⇀,ifp<0.5,D⇀′⋅ebl⋅cos2πl+X∗t,ifp≥0.5,where p is a random number in [0, 1].


### 2.3. Search for Prey

In nature, whales will go away from each other and randomly search for a new prey. In order to simulate the search behavior, a random search operation is used in the exploration phase when |A| ≥ 1. The current solution will update with a random solution in the current population instead of the best solution. The mathematical model is as follows:(7)D⇀=C⇀⋅X⇀randt−X⇀t,X⇀t+1=X⇀randt−A⇀⋅D⇀,where X⇀rand is a random solution in the current population.The pseudocode of the WOA is shown as Algorithm1.Algorithm 1: WOA.
Input: Population sizeNThe halting criterionOutput: The best solutionX∗Generate the initial populationXi (i = 1, 2, …, N)(1)
Calculate the fitness for each solution in the population(2)
Find the best solutionX∗ with the best fitness(3)
while the halting criterion is not satisfied do(4)
Updatea(5)
fori = 1 to Ndo(6)
UpdateA, C, l and p(7)
forj = 1 to mdo(8)
ifp<0.5then(9)
if (|A| < 1) then(10)
D=C⋅X∗j−Xij(11)
Xij=X∗j−A⋅D(12)
else if (|A| ≥ 1) then(13)
Select a random solutionXrand(14)
D=C⋅Xrandj−Xij(15)
Xij=Xrandj−A⋅D(16)
end if(17)
ifp≥0.5then(18)
D′=X∗j−Xij(19)
Xij=D′⋅ebl⋅cos2πl+X∗j(20)
end if(21)
end for(22)
Evaluate the fitness forXi(23)
end for(24)
Find the best solutionX∗(25)
end while


## 2.1. Encircling Prey

WOA supposes that the optimal solution with the best fitness in current population is the prey. All the other solutions will encircle the prey and update themselves toward the optimal solution. The mathematical model of encircling behavior is as follows:(1)D⇀=C⇀⋅X⇀∗t−X⇀t,(2)X⇀t+1=X⇀∗t−A⇀⋅D⇀,where t indicates the number of the current iteration, X is the current solution, X∗ is the best solution in the population, and A⇀ and C⇀ are the coefficient vectors calculated as follows:(3)A⇀=2a⇀⋅r⇀−a⇀,C⇀=2r⇀,where a⇀ is linearly decreased from 2 to 0 over the iterations and r⇀ is a random vector in [0, 1].


## 2.2. Bubble-Net Attacking Method

In order to simulate the behavior of bubble-net attacking, shrinking encircling mechanism and spiral updating position are designed as follows.


### 2.2.1. Shrinking Encircling Mechanism

The mathematical model of shrinking encircling mechanism is defined as equations (1) and (2). Because of a⇀ is linearly decreased from 2 to 0, the value of A⇀ is in the range of [−1, 1]. The current solution can update itself toward the best solution with a random step.


### 2.2.2. Spiral Updating Position

The new position of current solution updated with the best solution is obtained as follows:(4)D⇀′=X⇀∗t−X⇀t,(5)X⇀t+1=D⇀′⋅ebl⋅cos2πl+X∗t,where b is constant for defining logarithmic spiral and l is randomly generated number in the range [−1, 1]. To balance the two models, a selection parameter p is set. The selection parameter p in WOA is 0.5; there is a probability of 50% to choose either the shrinking encircling mechanism or the spiral updating position model. The finial mathematical model of bubble-net attacking is as follows:(6)X⇀t+1=X⇀∗t−A⇀⋅D⇀,ifp<0.5,D⇀′⋅ebl⋅cos2πl+X∗t,ifp≥0.5,where p is a random number in [0, 1].


## 2.2.1. Shrinking Encircling Mechanism

The mathematical model of shrinking encircling mechanism is defined as equations (1) and (2). Because of a⇀ is linearly decreased from 2 to 0, the value of A⇀ is in the range of [−1, 1]. The current solution can update itself toward the best solution with a random step.


## 2.2.2. Spiral Updating Position

The new position of current solution updated with the best solution is obtained as follows:(4)D⇀′=X⇀∗t−X⇀t,(5)X⇀t+1=D⇀′⋅ebl⋅cos2πl+X∗t,where b is constant for defining logarithmic spiral and l is randomly generated number in the range [−1, 1]. To balance the two models, a selection parameter p is set. The selection parameter p in WOA is 0.5; there is a probability of 50% to choose either the shrinking encircling mechanism or the spiral updating position model. The finial mathematical model of bubble-net attacking is as follows:(6)X⇀t+1=X⇀∗t−A⇀⋅D⇀,ifp<0.5,D⇀′⋅ebl⋅cos2πl+X∗t,ifp≥0.5,where p is a random number in [0, 1].


## 2.3. Search for Prey

In nature, whales will go away from each other and randomly search for a new prey. In order to simulate the search behavior, a random search operation is used in the exploration phase when |A| ≥ 1. The current solution will update with a random solution in the current population instead of the best solution. The mathematical model is as follows:(7)D⇀=C⇀⋅X⇀randt−X⇀t,X⇀t+1=X⇀randt−A⇀⋅D⇀,where X⇀rand is a random solution in the current population.The pseudocode of the WOA is shown as Algorithm1.Algorithm 1: WOA.
Input: Population sizeNThe halting criterionOutput: The best solutionX∗Generate the initial populationXi (i = 1, 2, …, N)(1)
Calculate the fitness for each solution in the population(2)
Find the best solutionX∗ with the best fitness(3)
while the halting criterion is not satisfied do(4)
Updatea(5)
fori = 1 to Ndo(6)
UpdateA, C, l and p(7)
forj = 1 to mdo(8)
ifp<0.5then(9)
if (|A| < 1) then(10)
D=C⋅X∗j−Xij(11)
Xij=X∗j−A⋅D(12)
else if (|A| ≥ 1) then(13)
Select a random solutionXrand(14)
D=C⋅Xrandj−Xij(15)
Xij=Xrandj−A⋅D(16)
end if(17)
ifp≥0.5then(18)
D′=X∗j−Xij(19)
Xij=D′⋅ebl⋅cos2πl+X∗j(20)
end if(21)
end for(22)
Evaluate the fitness forXi(23)
end for(24)
Find the best solutionX∗(25)
end while


## 3. The Proposed Approach

This section introduces a modified WOA named OCDWOA. Four strategies are introduced into the WOA. Opposition-based learning method (OBL) is used to initialize population and randomly update solutions in each iteration. Nonlinear parameter design is introduced to balance the ratio of exploitation phase and exploration phase. Density peak clustering strategy is adopted to decompose population into numbers of groups for enhancing the diversity of population. The effect of differential evolution (DE) is to exchange information between groups and improve the ability of jumping out of the local optimal solutions.


### 3.1. OCDWOA

The pseudocode of the OCDWOA algorithm is shown as Algorithm 2.Opposition-based learning is used to initialize the population (see lines 1–6 of Algorithm2) and to estimate the opposition solutions (see lines 39–46 of Algorithm 2). The fitness of each solution in the population is calculated and the best solution is found (see lines 7-8 in Algorithm 2). The parameters a, CR, A, C, l, and p are updated (see lines 10 and 13 in Algorithm 2). The population is decomposed into k groups by density peak clustering algorithm (see line 11 of Algorithm 2). The mutation operator of DE is introduced into the proposed algorithm (see lines 20–23 of Algorithm 2). A solution Xi is updated by r1 and r2 under mutation operator, where r1 is from the group including Xi and r2 is from the groups excluding Xi (see lines 14-15 of Algorithm 2). With the increase of iterations, the operation probability of DE declines. The basic three operators of WOA are used in lines 25-26, 29–31, and 34-35, respectively. The encircling prey is used in lines 25-26, the search for prey is introduced in lines 29–31, and the spiral updating position is applied in lines 34-35.Algorithm 2: OCDWOA.
Input: Population size NGroup sizekThe halting criterionOutput: The best solution X∗(1)
Generate the initial populationXi (i = 1, 2, …, N)(2)
fori = 1 to Ndo(3)
if (rand < 0.6)(4)
Generate new solution with OBL(5)
end if(6)
end for(7)
Calculate the fitness for each solution in the population(8)
Find the best solutionX∗ with the best fitness(9)
while the halting criterion is not satisfied do(10)
Updatea and CR(11)
Decompose population intok groups by density peak clustering algorithm(12)
fori = 1 to Ndo(13)
UpdateA, C, l and p(14)
R1 is the group with Xi and R2 are the groups without Xi(15)
Select randomlyr1 from R1 and r2 from R2(16)
jrand = randint (1, m)(17)
forj = 1 to mdo(18)
ifp<0.5then(19)
if (|A| < 1) then(20)
if rand (0, 1) ≤ CR andj == jrandthen(21)
UpdateF(22)
Select randomlyr1 from R1 and r2 from R2(23)
Uij=X∗j+FXr2j−Xr1j(24)
Else(25)
D=C⋅X∗j−Xij(26)
Xij=X∗j−A⋅D(27)
end if(28)
elseif (|A| ≥ 1) then(29)
Select a random solutionXrand(30)
D=C⋅Xrandj−Xij(31)
Xij=Xrandj−A⋅D(32)
end if(33)
ifp≥0.5then(34)
D′=X∗j−Xij(35)
Xij=D′⋅ebl⋅cos2πl+X∗j(36)
end if(37)
end if(38)
end for(39)
fori = 1 to Ndo(40)
if (rand < 0.6)(41)
Generate new solutionOi with OBL(42)
ifOi is better than Xithen(43)
Xi = Oi(44)
end if(45)
end if(46)
end for(47)
Find the best solutionX∗ with the best fitness(48)
end while


### 3.2. Opposition-Based Learning Method

It has been tested that opposition-based learning (OBL) method proposed by Rahnamayan et al. [35] is helpful for avoiding trapping into local optimal solutions.A solutionXi in d-dimension search space is denoted as Xi∗ = {Xi,1, Xi,2, Xi,3, …, Xi,d}. The opposition point Xi∗ of Xi can be calculated as(8)Xi,j′=Xub,j+Xlb,j−Xi,j,where Xub,j and Xlb,j denote the upper and lower bound of j-th dimension.OBL is used in the initialization step and solution update step at the end of each iteration. In the initialization step, OBL can generate a new solution on opposition side and enhance the diversity in the search space. In solution update step at the end of each iteration, OBL can generate a new solution on opposition side and choose the better solution between the new solution and the original solution. OBL can help WOA improve the search performance and jump out the local optimal solutions.


### 3.3. Nonlinear Parameter Design

An exponential functiona⇀ is introduced into the modified WOA. The value of coefficient vector A⇀ controls the proportion of exploitation phase and exploration phase. The nonlinear parameter design a⇀ is designed as follows [36]:(9)a=amin+amax−amin×e−2t/Max_iter,where amin and amax denote the minimum and maximum of a⇀, t is the number of the current iteration, and Max_iter is the maximum iteration.A large value ofa⇀ generates a large step size in exploitation phase, which can improve the searching speed. A small step size with a small value of a⇀ can be helpful for the exploration phase.


### 3.4. Density Peak Clustering

A density-based clustering strategy named density peak clustering algorithm [37] is used in the proposed WOA. The density peak clustering algorithm defines two important parameters for each solution i: the density ρi and distance δi. A scatter graph called decision diagram is constructed based on these two parameters. The decision diagram could decide the center of each clustering group and assign the remaining solutions to each clustering center.The mathematical model of the densityρi of the ith solution used in this paper is as follows:(10)ρi=∑je−dij/dc2,where j represents the j-th solution in current population, and i ≠ j, dij is the Euclidean distance between i-th and j-th solutions, and dc is a pre-set cut-off distance. A moderate value of dc could make the average number of neighbors around 1% to 2% of the population size.The distance parameterδi denotes the minimum distance among the i-th solution and the j-th solution, where the local density of j-th solution is higher than that of the i-th solution. The equation of δi is as follows:(11)δi=minj:ρj>ρidij.Figure1(a) shows the distribution of 28 solutions in two dimensions and Figure 1(b) is the decision diagram obtained by these 28 solutions. In Figure 1(b), the horizontal axis is ρ and the vertical axis is δ. The solutions with both large values of ρ and δ are defined as the clustering centers. In Figure 2, the 1st and 10th solutions with both large density and distance will be selected as the clustering centers. The 26th, 27th, and 28th solutions have large distance and low density, and the 7th and 8th solutions have large density and low distance. Therefore, these solutions cannot be selected as the clustering centers.Figure 1
(a) Solution distribution in two dimensions. (b) Decision diagram for each sample point.
(a)(b)Figure 2
Convergence graphs of OCDWOA with respect to other algorithms. (a)F1. (b) F2. (c) F3. (d) F4. (e) F5. (f) F6. (g) F7. (h) F8. (i) F9. (j) F10. (k) F11. (l) F12. (m) F13. (n) F14. (o) F15. (p) F16. (q) F17. (r) F18. (s) F19.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)(q)(r)(s)According to the rules of selecting the clustering centers, a modified equation of calculating quality of solutions is designed as follows:(12)Qi=ρi⋅δi.In this paper, a modified model of cut-off distancedc is introduced [38]. The transition distance Dis is calculated as follows:(13)Dis=1−e−dij2/2.The cut-off distancedc is calculated as follows:(14)dc=Dismax+Dismin2,where Dismax and Dismin are the maximum and minimum of transition distances in the current iteration.


### 3.5. Differential Evolution

For enhancing abilities of searching and jumping out of the local optimal solutions, the DE operator is introduced in this paper. The crossover and mutation operators in differential evolution (DE) [39] are calculated as follows:(15)Uij=X∗j+F×Xr2j−Xr1j,where Xr1 and Xr2 are the j-th dimension of the r1-th and r2-th solutions, respectively, U is the offspring, and F is the zoom factor. CR is the crossover probability updated as an exponential function:(16)CR=0.9×e−3t/Max_iter.


## 3.1. OCDWOA

The pseudocode of the OCDWOA algorithm is shown as Algorithm 2.Opposition-based learning is used to initialize the population (see lines 1–6 of Algorithm2) and to estimate the opposition solutions (see lines 39–46 of Algorithm 2). The fitness of each solution in the population is calculated and the best solution is found (see lines 7-8 in Algorithm 2). The parameters a, CR, A, C, l, and p are updated (see lines 10 and 13 in Algorithm 2). The population is decomposed into k groups by density peak clustering algorithm (see line 11 of Algorithm 2). The mutation operator of DE is introduced into the proposed algorithm (see lines 20–23 of Algorithm 2). A solution Xi is updated by r1 and r2 under mutation operator, where r1 is from the group including Xi and r2 is from the groups excluding Xi (see lines 14-15 of Algorithm 2). With the increase of iterations, the operation probability of DE declines. The basic three operators of WOA are used in lines 25-26, 29–31, and 34-35, respectively. The encircling prey is used in lines 25-26, the search for prey is introduced in lines 29–31, and the spiral updating position is applied in lines 34-35.Algorithm 2: OCDWOA.
Input: Population size NGroup sizekThe halting criterionOutput: The best solution X∗(1)
Generate the initial populationXi (i = 1, 2, …, N)(2)
fori = 1 to Ndo(3)
if (rand < 0.6)(4)
Generate new solution with OBL(5)
end if(6)
end for(7)
Calculate the fitness for each solution in the population(8)
Find the best solutionX∗ with the best fitness(9)
while the halting criterion is not satisfied do(10)
Updatea and CR(11)
Decompose population intok groups by density peak clustering algorithm(12)
fori = 1 to Ndo(13)
UpdateA, C, l and p(14)
R1 is the group with Xi and R2 are the groups without Xi(15)
Select randomlyr1 from R1 and r2 from R2(16)
jrand = randint (1, m)(17)
forj = 1 to mdo(18)
ifp<0.5then(19)
if (|A| < 1) then(20)
if rand (0, 1) ≤ CR andj == jrandthen(21)
UpdateF(22)
Select randomlyr1 from R1 and r2 from R2(23)
Uij=X∗j+FXr2j−Xr1j(24)
Else(25)
D=C⋅X∗j−Xij(26)
Xij=X∗j−A⋅D(27)
end if(28)
elseif (|A| ≥ 1) then(29)
Select a random solutionXrand(30)
D=C⋅Xrandj−Xij(31)
Xij=Xrandj−A⋅D(32)
end if(33)
ifp≥0.5then(34)
D′=X∗j−Xij(35)
Xij=D′⋅ebl⋅cos2πl+X∗j(36)
end if(37)
end if(38)
end for(39)
fori = 1 to Ndo(40)
if (rand < 0.6)(41)
Generate new solutionOi with OBL(42)
ifOi is better than Xithen(43)
Xi = Oi(44)
end if(45)
end if(46)
end for(47)
Find the best solutionX∗ with the best fitness(48)
end while


## 3.2. Opposition-Based Learning Method

It has been tested that opposition-based learning (OBL) method proposed by Rahnamayan et al. [35] is helpful for avoiding trapping into local optimal solutions.A solutionXi in d-dimension search space is denoted as Xi∗ = {Xi,1, Xi,2, Xi,3, …, Xi,d}. The opposition point Xi∗ of Xi can be calculated as(8)Xi,j′=Xub,j+Xlb,j−Xi,j,where Xub,j and Xlb,j denote the upper and lower bound of j-th dimension.OBL is used in the initialization step and solution update step at the end of each iteration. In the initialization step, OBL can generate a new solution on opposition side and enhance the diversity in the search space. In solution update step at the end of each iteration, OBL can generate a new solution on opposition side and choose the better solution between the new solution and the original solution. OBL can help WOA improve the search performance and jump out the local optimal solutions.


## 3.3. Nonlinear Parameter Design

An exponential functiona⇀ is introduced into the modified WOA. The value of coefficient vector A⇀ controls the proportion of exploitation phase and exploration phase. The nonlinear parameter design a⇀ is designed as follows [36]:(9)a=amin+amax−amin×e−2t/Max_iter,where amin and amax denote the minimum and maximum of a⇀, t is the number of the current iteration, and Max_iter is the maximum iteration.A large value ofa⇀ generates a large step size in exploitation phase, which can improve the searching speed. A small step size with a small value of a⇀ can be helpful for the exploration phase.


## 3.4. Density Peak Clustering

A density-based clustering strategy named density peak clustering algorithm [37] is used in the proposed WOA. The density peak clustering algorithm defines two important parameters for each solution i: the density ρi and distance δi. A scatter graph called decision diagram is constructed based on these two parameters. The decision diagram could decide the center of each clustering group and assign the remaining solutions to each clustering center.The mathematical model of the densityρi of the ith solution used in this paper is as follows:(10)ρi=∑je−dij/dc2,where j represents the j-th solution in current population, and i ≠ j, dij is the Euclidean distance between i-th and j-th solutions, and dc is a pre-set cut-off distance. A moderate value of dc could make the average number of neighbors around 1% to 2% of the population size.The distance parameterδi denotes the minimum distance among the i-th solution and the j-th solution, where the local density of j-th solution is higher than that of the i-th solution. The equation of δi is as follows:(11)δi=minj:ρj>ρidij.Figure1(a) shows the distribution of 28 solutions in two dimensions and Figure 1(b) is the decision diagram obtained by these 28 solutions. In Figure 1(b), the horizontal axis is ρ and the vertical axis is δ. The solutions with both large values of ρ and δ are defined as the clustering centers. In Figure 2, the 1st and 10th solutions with both large density and distance will be selected as the clustering centers. The 26th, 27th, and 28th solutions have large distance and low density, and the 7th and 8th solutions have large density and low distance. Therefore, these solutions cannot be selected as the clustering centers.Figure 1
(a) Solution distribution in two dimensions. (b) Decision diagram for each sample point.
(a)(b)Figure 2
Convergence graphs of OCDWOA with respect to other algorithms. (a)F1. (b) F2. (c) F3. (d) F4. (e) F5. (f) F6. (g) F7. (h) F8. (i) F9. (j) F10. (k) F11. (l) F12. (m) F13. (n) F14. (o) F15. (p) F16. (q) F17. (r) F18. (s) F19.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)(q)(r)(s)According to the rules of selecting the clustering centers, a modified equation of calculating quality of solutions is designed as follows:(12)Qi=ρi⋅δi.In this paper, a modified model of cut-off distancedc is introduced [38]. The transition distance Dis is calculated as follows:(13)Dis=1−e−dij2/2.The cut-off distancedc is calculated as follows:(14)dc=Dismax+Dismin2,where Dismax and Dismin are the maximum and minimum of transition distances in the current iteration.


## 3.5. Differential Evolution

For enhancing abilities of searching and jumping out of the local optimal solutions, the DE operator is introduced in this paper. The crossover and mutation operators in differential evolution (DE) [39] are calculated as follows:(15)Uij=X∗j+F×Xr2j−Xr1j,where Xr1 and Xr2 are the j-th dimension of the r1-th and r2-th solutions, respectively, U is the offspring, and F is the zoom factor. CR is the crossover probability updated as an exponential function:(16)CR=0.9×e−3t/Max_iter.


## 4. Experiment Results and Discussion

### 4.1. Parameter Settings and Benchmark Functions

In this section, in order to analyze the performance of the proposed algorithm, 19 benchmark functions [40] F1–F19 have been used. The expressions of all the functions are given in Table 1. F1–F7 are the unimodal functions with only one global optimal solution. F8–F12 are the multimodal functions with many local optimal solutions. These functions can test the performance of jumping out the local optimal solutions. F13–F19 are the fixed-dimension multimodal functions, which can evaluate the convergence property of an algorithm.Table 1
Description of benchmark functions.
FunctionDimRangefminUnimodal functionsF1=∑i=1nxi230[−100, 100]0F2=∑i=1nxi+∏i=1nxi30[−10, 10]0F3=∑i=1n∑j=1ixj230[−100, 100]0F4=maxixi,1≤i≤n30[−100, 100]0F5=∑i=1n−1100xi+1−xi22+xi−1230[−30, 30]0F6=∑i=1nxi+0.5230[−100, 100]0F7=∑i=1nixi4+random0,130[−1.28, 1.28]0Multimodal functionsF8=∑i=1nxi2−10cos2πxi+1030[−5.12.5.12]0F9=−20exp−0.21/n∑i=1nxi2−exp1/n∑i=1ncos2πxi+20+e30[−32, 32]0F10=1/4000∑i=1nxi2+∏i=1ncosxi/i+130[−600, 600]0F11=π/n10sinπy1+∑i=1n−1yi−121+10sin2πyi+1+yn−12+∑i=1nuxi,10,100,4yi=1+xi+1/4uxi,a,k,m=kxi−ama<xi0−a<xi<ak−xi−amxi<−a30[−50, 50]0F12=0.1sin23πx1+∑i=1n−1xi−121+sin23πxi+1+xn−121+sin22πxn+∑i=1nuxi,5,100,430[−50, 50]0Fixed-dimension multimodal functionsF13=1/500+∑j=1251/j+∑i=12xi−aij6−12[−65, 65]0F14=∑i=111ai+xibi2+bix2/bi2+bix3+x424[−5, 5]0.00030F15=4x12−2.1x14+1/3x16+x1x2−4x22+4x242[−5, 5]−1.0316F16=x2−5.1/4π2x12+5/πx1−62+101−1/8πcosx1+102[−5, 5]0.398F17=1+x1+x2+1219−14x1+3x12−14x2+6x1x2+3x22×30+2x1−3x22×18−32x1+12x12+48x2−36x1x2+27x222[−2, 2]3F18=−∑i=14ciexp−∑j=13aijxj−pij23[1, 3]−3.86F19=−∑i=14ciexp−∑j=16aijxj−pij26[0, 3]−3.12The proposed OCDWOA is compared with the classical WOA and its typical variants, such as OWOA [41], OEWOA [42], and IWOA [43]. The parameter b is set to 1 for all the algorithms. The parameter a⇀ decreases linearly from 2 to 0 for WOA, OWOA, and IWOA. For OEWOA, the parameter a⇀ nonlinearly decreases from 2 to 0. For IWOA, the crossover probability CR is set to 0.9, and the zoom factor F is random within (0.2, 0.8).For all the algorithms, population size is 30 and maximum iteration is 500. The experiment is repeated 30 times for each benchmark function.


### 4.2. Experiment Results

The statistics of the experimental results is demonstrated in Table2. The best fitness, mean fitness, and standard deviation of OCDWOA and each compared algorithm are given.Table 2
Results of OCDWOA and compared algorithms.
FunctionAlgorithmBestMeanStdF1WOA3.45E − 841.84E − 746.85E − 74OWOA3.67E − 871.54E − 726.15E − 72OEWOA8.10E − 1043.41E − 871.61E − 86IWOA3.88E − 314.18E − 281.52E − 27OCDWOA1.16E − 942.68E − 828.18E − 82F2WOA3.76E − 594.29E − 492.30E − 48OWOA8.72E − 583.71E − 509.41E − 50OEWOA8.88E − 682.26E − 591.15E − 58IWOA6.02E − 206.36E − 196.72E − 19OCDWOA1.85E − 631.86E − 564.92E − 56F3WOA29331.3904346962.3122811654.1269OWOA23102.6358545115.12149775.63579OEWOA10896.0878737765.4224912462.65209IWOA0.02049194.76100185112.20526468OCDWOA9123.14658930378.8800311250.945F4WOA2.84594453553.9849200928.64604813OWOA0.22168328354.5971715128.7671405OEWOA0.39550066651.0564977130.09654391IWOA2.53E − 051.43E − 031.46E − 03OCDWOA0.09125994227.1570554822.72468237F5WOA27.0509657927.982064580.492345074OWOA27.046271327.86281190.3574981OEWOA27.4796158428.172559650.377986162IWOA24.7260899126.219984581.236007685OCDWOA26.9433084927.752361940.519791973F6WOA0.0393524230.3081862760.178197178OWOA0.0589002570.2294445010.138575416OEWOA0.0942035930.4719325210.211010397IWOA0.1325538970.6348501880.363709611OCDWOA0.0118807430.1777461210.202309313F7WOA7.31E − 053.49E − 034.08E − 03OWOA9.68E − 061.20E − 032.15E − 03OEWOA1.78E − 058.51E − 049.84E − 04IWOA0.0011707660.0045713040.003016255OCDWOA1.19E − 058.98E − 041.21E − 03F8WOA05.68E − 152.29E − 14OWOA000OEWOA000IWOA022.3668107124.91175078OCDWOA000F9WOA8.88E − 164.32E − 152.55E − 15OWOA8.88E − 163.85E − 152.48E − 15OEWOA8.88E − 164.20E − 152.27E − 15IWOA2.22E − 144.29E − 141.27E − 14OCDWOA8.88E − 164.44E − 152.09E − 15F10WOA00.0142391430.055902557OWOA00.0063693050.034886118OEWOA000IWOA00.0056015190.009188179OCDWOA00.003071750.016824667F11WOA0.0033033070.0207038830.013561528OWOA0.0028863910.0115234240.008134223OEWOA0.0088326440.0229317860.008957016IWOA0.012898710.041217510.023964971OCDWOA0.0009847440.0080862270.00817336F12WOA0.0917596920.4769726390.2731912OWOA0.0712933810.3439783840.176667484OEWOA0.2645983220.5587882510.219290838IWOA0.3600523840.8261869960.250307701OCDWOA0.0089727270.3150169920.165937507F13WOA0.9980038383.0877078173.643013869OWOA0.9980038381.5274928330.770156088OEWOA0.9980038381.9910646761.244254244IWOA0.9980038385.4916273484.859439425OCDWOA0.9980038381.3288191510.7056337F14WOA0.0003093540.0007547010.000591433OWOA0.0003080460.0006927850.000506369OEWOA0.0003080420.0005994120.000433827IWOA0.0003074940.0057638640.008959844OCDWOA0.0003080150.000683180.000396378F15WOA−1.031628453−1.0316284537.57E − 10OWOA−1.031628453−1.0316284531.10E − 09OEWOA−1.031628453−1.0316284535.29E − 10IWOA−1.031628453−1.0316284532.29E − 10OCDWOA−1.031628453−1.0316284533.68E − 11F16WOA0.3978873580.3978944121.60E − 05OWOA0.3978873580.3978997212.59E − 05OEWOA0.3978873580.3979169920.000107977IWOA0.3978873580.3978875855.13E − 07OCDWOA0.3978873580.3978885162.38E − 06F17WOA3.0000000013.0000438727.74E − 05OWOA3.000000023.0000446350.000152146OEWOA3.0000000013.0001236650.000226823IWOA38.40000004220.55035867OCDWOA3.000000013.0000063551.60E − 05F18WOA−3.862688519−3.8572620650.007143618OWOA−3.862781863−3.8618154230.003572799OEWOA−3.862782142−3.857266970.019011496IWOA−3.862782147−3.8625215550.001426916OCDWOA−3.862782147−3.8627759169.47E − 06F19WOA−3.321580972−3.2235978770.10220952OWOA−3.321545137−3.2441433490.085289027OEWOA−3.320833167−3.2130348030.091087411IWOA−3.32199476−3.2700869440.060385879OCDWOA−3.321994702−3.2783247840.058345799The bold values mean the best results.ForF6, F9, F11, F12, F13, F15, F17, and F18, the performance of OCDWOA is the absolutely best. For F14, the results obtained by all the variants of WOA are comparable. For F1, F2, F3, F4, F5, F10, and F16, OCDWOA obtains the second rank. For the simple functions F1 and F2 without any local optimal solution, OCDWOA is worse than OEWOA. The DE operator used in OCDWOA can enhance the search performance. The convergence speed would be reduced, though the operation probability of DE declines in exploration phase. For the nonseparable unimodal functions F3, F4, and F5, OCDWOA is only worse than IWOA. For the multimodal function F16, the result reveals the positive influence of DE.ForF7 and F8, the performance of OCDWOA is similar to ones of OWOA and OEWOA. For F10, OCDWOA gets the second rank. For F8 and F10, IWOA can get the optimal solution. However, IWOA fails in solving the function F7. The mean value and standard deviation value of IWOA are also not ideal. It shows that the DE operator is helpful for jumping out the local optimal solutions, but the stability may be impaired for the multimodal functions with easy or few local optimal solutions. In OCDWOA, the clustering operator and reducing the operation probability of DE can drop off the influence of DE. The DE operator can improve the performance of jump out the local optimal solutions. The DE operator and OBL can enhance the search performance of algorithm. Dynamical adjustment of a⇀ can avoid premature convergence. The clustering operator enhances information exchange between solutions.The results analyzed by the Wilcoxon rank-sum test are shown in Table3. In the Wilcoxon rank-sum test, a significance level is set to 0.05. In Table 3, “+,” “−,” and “ = ” indicate that the result of the proposed algorithm OCDWOA is significantly better, significantly worse, and statistically similar to that of the compared algorithm, respectively. For most test functions, the performance of OCDWOA is significantly better or statistically similar to the other variants of WOA.Table 3
The Wilcoxon rank-sum test results of OCDWOA and compared algorithms.
OCDWOA-WOAOCDWOA-OWOAOCDWOA-OEWOAOCDWOA-IWOAComparison+/=/−+/=/−+/=/−+/=/−Result14/5/010/9/011/6/210/4/5The convergence curves and boxplots are shown in Figures2 and 3, respectively.Figure 3
Boxplots of OCDWOA with respect to other algorithms. (a)F1. (b) F2. (c) F3. (d) F4. (e) F5. (f) F6. (g) F7. (h) F8. (i) F9. (j) F10. (k) F11. (l) F12. (m) F13. (n) F15. (o) F16. (p) F17. (q) F18. (r) F19. (s) F19.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)(q)(r)(s)In Figure2, for most functions, the convergence curves of OCDWOA are better or statistically similar to ones of the other variants of WOA. For F3 and F4, the convergence curves of OCDWOA are stabler than those of IWOA, but the convergence rates of OCDWOA are slower than those of IWOA. In OCDWOA, the solutions selected in the crossover and mutation operators of DE are based on the clustering. The clustering operator reduces the randomness of the DE operator. The stability of OCDWOA has been improved.In Figure3, for most functions, the stability of OCDWOA is better or statistically similar to the stability of the other variants of WOA. Only for F3 and F4, the stability of OCDWOA is worse than the stability of IWOA.


### 4.3. Application of Seismic Inversion Problems

#### 4.3.1. Seismic Inversion Problem

Seismic inversion is the key technology of quantitative reservoir prediction and prestack waveform inversion [44] is a newly developed seismic inversion technique. In prestack waveform inversion technique, the complex wave propagation effect is considered and the complete waveform information is fully used. The prestack waveform inversion method assumes that each CMP/CDP point of underground medium layer structure with local one-dimensional or level uses the reflectivity method forward prestack data [45]. The prestack waveform inversion method searches the best fitting with the actual prestack seismic data of underground through genetic algorithm, simulated annealing, and other global optimization algorithms [46]. At present, there are many problems in prestack waveform inversion technology that need to be further studied and the resolution of prestack waveform inversion technology needs to be improved. The existing optimization algorithms used in prestack waveform inversion are mainly traditional intelligent optimization algorithms. These algorithms are prone to premature convergence and local optimal solution when dealing with large data volume, highly nonlinear, multiparameter, and multiextremum prestack waveform inversion problems.In this paper, a theoretical geological model with 11 layers is established [47]. Each layer in the model contains 3 parameters: velocity of longitudinal wave (VP), velocity of shear wave (VS), and density (ρ). The values of model are shown in Table 4.Table 4
Values of the theoretical geological model.
VP (m·s−1)VS (m·s−1)ρ (g·cm−3)Layer 119808512Layer 2218516052.2Layer 319838562Layer 4220016222.2Layer 5243715102.3Layer 6266313802.4Layer 719838512Layer 8243715012.3Layer 919808562Layer 10265613802.4Layer 1119838512In this paper, the reflection coefficients of the model will be calculated by Zoeppritz equation as follows:(17)RPRSTPTS=−sinθ1−cosϕ1sinθ2cosϕ2cosθ1−sinϕ1cosθ2−sinϕ2sin2θ1VP1VS1cos2ϕ1ρ2VS22VP1ρ1VS12VP2cos2ϕ1ρ2VS2VP1ρ1VS12cos2ϕ2−cos2ϕ1VS1VP1sin2ϕ1ρ2VP2ρ1VP1cos2ϕ2−ρ2VS2ρ1VP1sin2ϕ2−1sinθ1cosθ1sin2θ1cos2ϕ1,where VP1, VS1, and ρ1 are the parameters of the up-layer, VP2, VS2, and ρ2 are the parameters of the low-layer, θ1 and θ2 are the incident angle and the refraction angle of the P-wave, respectively, and Φ1 and Φ2 are the incident angle and the refraction angle of the S-wave, respectively.Φ1, θ2, and Φ2 are calculated as follows, respectively:(18)ϕ1=arcsinθ1×VS1VP1,θ2=arcsinθ1×VS2VP1,ϕ2=arcsinθ1×VP2VP1.Incident anglesθ1 of 50 reflection coefficients are from 1° to 50°. The 50 seismic records are obtained by convolving the reflection coefficients with the Ricker wavelet with 35 Hz as follows:(19)S=R∗Rw,where S is the seismic record, R is the reflection coefficient, and Rw is the ricker wavelet (35 Hz).The 50 seismic records will be divided into three groups. The three groups are small incident angle group (Group 1), medium incident angle group (Group 2), and large incident angle group (Group 3). Group 1 consists of the seismic records whose incident angles are from 1° to 15°. Group 2 consists of the seismic records whose incident angles are from 16° to 32°. Group 3 consists of the seismic records whose incident angles are from 33° to 50°. Superimposed seismic records S1, S2, and S3 are obtained by synthesizing three groups of seismic records, respectively. The superimposed seismic records S1(X), S2(X), and S3(X) of the solution X are obtained the same way as the model.The fitness of the solutionX is expressed as follows:(20)FSIX=13S1X−S12+S2X−S22+S3X−S32.


#### 4.3.2. Experiment Results

The experiment is repeated 10 times for each algorithm. Table5 shows the best fitness, mean fitness, and standard deviation of test results on seismic inversion problem. The convergence curves and boxplots are shown in Figure 4. From the test results, the performance of OCDWOA for seismic inversion problem is the best.Table 5
Results of different algorithms.
ProblemAlgorithmBestMeanStdSeismic inversionWOA0.0071787540.01241271160.0038698430OWOA0.0061682620.01195090700.0053173170OEWOA0.0069015790.01120280110.0035417924IWOA0.0061365810.01289760650.0041582246OCDWOA0.0026622290.00559751190.0019518514The bold values mean the best results.Figure 4
Results of OCDWOA with respect to other algorithms. (a) Convergence graph of the algorithms. (b) Boxplots of the algorithms.
(a)(b)In Table5, the best fitness, mean fitness, and standard deviation of OCDWOA are all obviously smaller than those of the comparative algorithms. In Figure 4(a), the convergence rate of OCDWOA is better than those of the comparative algorithms. The convergence curve of OCDWOA is stable. In Figure 4(b), the stability of OCDWOA is the best among the comparative algorithms.


## 4.1. Parameter Settings and Benchmark Functions

In this section, in order to analyze the performance of the proposed algorithm, 19 benchmark functions [40] F1–F19 have been used. The expressions of all the functions are given in Table 1. F1–F7 are the unimodal functions with only one global optimal solution. F8–F12 are the multimodal functions with many local optimal solutions. These functions can test the performance of jumping out the local optimal solutions. F13–F19 are the fixed-dimension multimodal functions, which can evaluate the convergence property of an algorithm.Table 1
Description of benchmark functions.
FunctionDimRangefminUnimodal functionsF1=∑i=1nxi230[−100, 100]0F2=∑i=1nxi+∏i=1nxi30[−10, 10]0F3=∑i=1n∑j=1ixj230[−100, 100]0F4=maxixi,1≤i≤n30[−100, 100]0F5=∑i=1n−1100xi+1−xi22+xi−1230[−30, 30]0F6=∑i=1nxi+0.5230[−100, 100]0F7=∑i=1nixi4+random0,130[−1.28, 1.28]0Multimodal functionsF8=∑i=1nxi2−10cos2πxi+1030[−5.12.5.12]0F9=−20exp−0.21/n∑i=1nxi2−exp1/n∑i=1ncos2πxi+20+e30[−32, 32]0F10=1/4000∑i=1nxi2+∏i=1ncosxi/i+130[−600, 600]0F11=π/n10sinπy1+∑i=1n−1yi−121+10sin2πyi+1+yn−12+∑i=1nuxi,10,100,4yi=1+xi+1/4uxi,a,k,m=kxi−ama<xi0−a<xi<ak−xi−amxi<−a30[−50, 50]0F12=0.1sin23πx1+∑i=1n−1xi−121+sin23πxi+1+xn−121+sin22πxn+∑i=1nuxi,5,100,430[−50, 50]0Fixed-dimension multimodal functionsF13=1/500+∑j=1251/j+∑i=12xi−aij6−12[−65, 65]0F14=∑i=111ai+xibi2+bix2/bi2+bix3+x424[−5, 5]0.00030F15=4x12−2.1x14+1/3x16+x1x2−4x22+4x242[−5, 5]−1.0316F16=x2−5.1/4π2x12+5/πx1−62+101−1/8πcosx1+102[−5, 5]0.398F17=1+x1+x2+1219−14x1+3x12−14x2+6x1x2+3x22×30+2x1−3x22×18−32x1+12x12+48x2−36x1x2+27x222[−2, 2]3F18=−∑i=14ciexp−∑j=13aijxj−pij23[1, 3]−3.86F19=−∑i=14ciexp−∑j=16aijxj−pij26[0, 3]−3.12The proposed OCDWOA is compared with the classical WOA and its typical variants, such as OWOA [41], OEWOA [42], and IWOA [43]. The parameter b is set to 1 for all the algorithms. The parameter a⇀ decreases linearly from 2 to 0 for WOA, OWOA, and IWOA. For OEWOA, the parameter a⇀ nonlinearly decreases from 2 to 0. For IWOA, the crossover probability CR is set to 0.9, and the zoom factor F is random within (0.2, 0.8).For all the algorithms, population size is 30 and maximum iteration is 500. The experiment is repeated 30 times for each benchmark function.


## 4.2. Experiment Results

The statistics of the experimental results is demonstrated in Table2. The best fitness, mean fitness, and standard deviation of OCDWOA and each compared algorithm are given.Table 2
Results of OCDWOA and compared algorithms.
FunctionAlgorithmBestMeanStdF1WOA3.45E − 841.84E − 746.85E − 74OWOA3.67E − 871.54E − 726.15E − 72OEWOA8.10E − 1043.41E − 871.61E − 86IWOA3.88E − 314.18E − 281.52E − 27OCDWOA1.16E − 942.68E − 828.18E − 82F2WOA3.76E − 594.29E − 492.30E − 48OWOA8.72E − 583.71E − 509.41E − 50OEWOA8.88E − 682.26E − 591.15E − 58IWOA6.02E − 206.36E − 196.72E − 19OCDWOA1.85E − 631.86E − 564.92E − 56F3WOA29331.3904346962.3122811654.1269OWOA23102.6358545115.12149775.63579OEWOA10896.0878737765.4224912462.65209IWOA0.02049194.76100185112.20526468OCDWOA9123.14658930378.8800311250.945F4WOA2.84594453553.9849200928.64604813OWOA0.22168328354.5971715128.7671405OEWOA0.39550066651.0564977130.09654391IWOA2.53E − 051.43E − 031.46E − 03OCDWOA0.09125994227.1570554822.72468237F5WOA27.0509657927.982064580.492345074OWOA27.046271327.86281190.3574981OEWOA27.4796158428.172559650.377986162IWOA24.7260899126.219984581.236007685OCDWOA26.9433084927.752361940.519791973F6WOA0.0393524230.3081862760.178197178OWOA0.0589002570.2294445010.138575416OEWOA0.0942035930.4719325210.211010397IWOA0.1325538970.6348501880.363709611OCDWOA0.0118807430.1777461210.202309313F7WOA7.31E − 053.49E − 034.08E − 03OWOA9.68E − 061.20E − 032.15E − 03OEWOA1.78E − 058.51E − 049.84E − 04IWOA0.0011707660.0045713040.003016255OCDWOA1.19E − 058.98E − 041.21E − 03F8WOA05.68E − 152.29E − 14OWOA000OEWOA000IWOA022.3668107124.91175078OCDWOA000F9WOA8.88E − 164.32E − 152.55E − 15OWOA8.88E − 163.85E − 152.48E − 15OEWOA8.88E − 164.20E − 152.27E − 15IWOA2.22E − 144.29E − 141.27E − 14OCDWOA8.88E − 164.44E − 152.09E − 15F10WOA00.0142391430.055902557OWOA00.0063693050.034886118OEWOA000IWOA00.0056015190.009188179OCDWOA00.003071750.016824667F11WOA0.0033033070.0207038830.013561528OWOA0.0028863910.0115234240.008134223OEWOA0.0088326440.0229317860.008957016IWOA0.012898710.041217510.023964971OCDWOA0.0009847440.0080862270.00817336F12WOA0.0917596920.4769726390.2731912OWOA0.0712933810.3439783840.176667484OEWOA0.2645983220.5587882510.219290838IWOA0.3600523840.8261869960.250307701OCDWOA0.0089727270.3150169920.165937507F13WOA0.9980038383.0877078173.643013869OWOA0.9980038381.5274928330.770156088OEWOA0.9980038381.9910646761.244254244IWOA0.9980038385.4916273484.859439425OCDWOA0.9980038381.3288191510.7056337F14WOA0.0003093540.0007547010.000591433OWOA0.0003080460.0006927850.000506369OEWOA0.0003080420.0005994120.000433827IWOA0.0003074940.0057638640.008959844OCDWOA0.0003080150.000683180.000396378F15WOA−1.031628453−1.0316284537.57E − 10OWOA−1.031628453−1.0316284531.10E − 09OEWOA−1.031628453−1.0316284535.29E − 10IWOA−1.031628453−1.0316284532.29E − 10OCDWOA−1.031628453−1.0316284533.68E − 11F16WOA0.3978873580.3978944121.60E − 05OWOA0.3978873580.3978997212.59E − 05OEWOA0.3978873580.3979169920.000107977IWOA0.3978873580.3978875855.13E − 07OCDWOA0.3978873580.3978885162.38E − 06F17WOA3.0000000013.0000438727.74E − 05OWOA3.000000023.0000446350.000152146OEWOA3.0000000013.0001236650.000226823IWOA38.40000004220.55035867OCDWOA3.000000013.0000063551.60E − 05F18WOA−3.862688519−3.8572620650.007143618OWOA−3.862781863−3.8618154230.003572799OEWOA−3.862782142−3.857266970.019011496IWOA−3.862782147−3.8625215550.001426916OCDWOA−3.862782147−3.8627759169.47E − 06F19WOA−3.321580972−3.2235978770.10220952OWOA−3.321545137−3.2441433490.085289027OEWOA−3.320833167−3.2130348030.091087411IWOA−3.32199476−3.2700869440.060385879OCDWOA−3.321994702−3.2783247840.058345799The bold values mean the best results.ForF6, F9, F11, F12, F13, F15, F17, and F18, the performance of OCDWOA is the absolutely best. For F14, the results obtained by all the variants of WOA are comparable. For F1, F2, F3, F4, F5, F10, and F16, OCDWOA obtains the second rank. For the simple functions F1 and F2 without any local optimal solution, OCDWOA is worse than OEWOA. The DE operator used in OCDWOA can enhance the search performance. The convergence speed would be reduced, though the operation probability of DE declines in exploration phase. For the nonseparable unimodal functions F3, F4, and F5, OCDWOA is only worse than IWOA. For the multimodal function F16, the result reveals the positive influence of DE.ForF7 and F8, the performance of OCDWOA is similar to ones of OWOA and OEWOA. For F10, OCDWOA gets the second rank. For F8 and F10, IWOA can get the optimal solution. However, IWOA fails in solving the function F7. The mean value and standard deviation value of IWOA are also not ideal. It shows that the DE operator is helpful for jumping out the local optimal solutions, but the stability may be impaired for the multimodal functions with easy or few local optimal solutions. In OCDWOA, the clustering operator and reducing the operation probability of DE can drop off the influence of DE. The DE operator can improve the performance of jump out the local optimal solutions. The DE operator and OBL can enhance the search performance of algorithm. Dynamical adjustment of a⇀ can avoid premature convergence. The clustering operator enhances information exchange between solutions.The results analyzed by the Wilcoxon rank-sum test are shown in Table3. In the Wilcoxon rank-sum test, a significance level is set to 0.05. In Table 3, “+,” “−,” and “ = ” indicate that the result of the proposed algorithm OCDWOA is significantly better, significantly worse, and statistically similar to that of the compared algorithm, respectively. For most test functions, the performance of OCDWOA is significantly better or statistically similar to the other variants of WOA.Table 3
The Wilcoxon rank-sum test results of OCDWOA and compared algorithms.
OCDWOA-WOAOCDWOA-OWOAOCDWOA-OEWOAOCDWOA-IWOAComparison+/=/−+/=/−+/=/−+/=/−Result14/5/010/9/011/6/210/4/5The convergence curves and boxplots are shown in Figures2 and 3, respectively.Figure 3
Boxplots of OCDWOA with respect to other algorithms. (a)F1. (b) F2. (c) F3. (d) F4. (e) F5. (f) F6. (g) F7. (h) F8. (i) F9. (j) F10. (k) F11. (l) F12. (m) F13. (n) F15. (o) F16. (p) F17. (q) F18. (r) F19. (s) F19.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)(q)(r)(s)In Figure2, for most functions, the convergence curves of OCDWOA are better or statistically similar to ones of the other variants of WOA. For F3 and F4, the convergence curves of OCDWOA are stabler than those of IWOA, but the convergence rates of OCDWOA are slower than those of IWOA. In OCDWOA, the solutions selected in the crossover and mutation operators of DE are based on the clustering. The clustering operator reduces the randomness of the DE operator. The stability of OCDWOA has been improved.In Figure3, for most functions, the stability of OCDWOA is better or statistically similar to the stability of the other variants of WOA. Only for F3 and F4, the stability of OCDWOA is worse than the stability of IWOA.


## 4.3. Application of Seismic Inversion Problems

### 4.3.1. Seismic Inversion Problem

Seismic inversion is the key technology of quantitative reservoir prediction and prestack waveform inversion [44] is a newly developed seismic inversion technique. In prestack waveform inversion technique, the complex wave propagation effect is considered and the complete waveform information is fully used. The prestack waveform inversion method assumes that each CMP/CDP point of underground medium layer structure with local one-dimensional or level uses the reflectivity method forward prestack data [45]. The prestack waveform inversion method searches the best fitting with the actual prestack seismic data of underground through genetic algorithm, simulated annealing, and other global optimization algorithms [46]. At present, there are many problems in prestack waveform inversion technology that need to be further studied and the resolution of prestack waveform inversion technology needs to be improved. The existing optimization algorithms used in prestack waveform inversion are mainly traditional intelligent optimization algorithms. These algorithms are prone to premature convergence and local optimal solution when dealing with large data volume, highly nonlinear, multiparameter, and multiextremum prestack waveform inversion problems.In this paper, a theoretical geological model with 11 layers is established [47]. Each layer in the model contains 3 parameters: velocity of longitudinal wave (VP), velocity of shear wave (VS), and density (ρ). The values of model are shown in Table 4.Table 4
Values of the theoretical geological model.
VP (m·s−1)VS (m·s−1)ρ (g·cm−3)Layer 119808512Layer 2218516052.2Layer 319838562Layer 4220016222.2Layer 5243715102.3Layer 6266313802.4Layer 719838512Layer 8243715012.3Layer 919808562Layer 10265613802.4Layer 1119838512In this paper, the reflection coefficients of the model will be calculated by Zoeppritz equation as follows:(17)RPRSTPTS=−sinθ1−cosϕ1sinθ2cosϕ2cosθ1−sinϕ1cosθ2−sinϕ2sin2θ1VP1VS1cos2ϕ1ρ2VS22VP1ρ1VS12VP2cos2ϕ1ρ2VS2VP1ρ1VS12cos2ϕ2−cos2ϕ1VS1VP1sin2ϕ1ρ2VP2ρ1VP1cos2ϕ2−ρ2VS2ρ1VP1sin2ϕ2−1sinθ1cosθ1sin2θ1cos2ϕ1,where VP1, VS1, and ρ1 are the parameters of the up-layer, VP2, VS2, and ρ2 are the parameters of the low-layer, θ1 and θ2 are the incident angle and the refraction angle of the P-wave, respectively, and Φ1 and Φ2 are the incident angle and the refraction angle of the S-wave, respectively.Φ1, θ2, and Φ2 are calculated as follows, respectively:(18)ϕ1=arcsinθ1×VS1VP1,θ2=arcsinθ1×VS2VP1,ϕ2=arcsinθ1×VP2VP1.Incident anglesθ1 of 50 reflection coefficients are from 1° to 50°. The 50 seismic records are obtained by convolving the reflection coefficients with the Ricker wavelet with 35 Hz as follows:(19)S=R∗Rw,where S is the seismic record, R is the reflection coefficient, and Rw is the ricker wavelet (35 Hz).The 50 seismic records will be divided into three groups. The three groups are small incident angle group (Group 1), medium incident angle group (Group 2), and large incident angle group (Group 3). Group 1 consists of the seismic records whose incident angles are from 1° to 15°. Group 2 consists of the seismic records whose incident angles are from 16° to 32°. Group 3 consists of the seismic records whose incident angles are from 33° to 50°. Superimposed seismic records S1, S2, and S3 are obtained by synthesizing three groups of seismic records, respectively. The superimposed seismic records S1(X), S2(X), and S3(X) of the solution X are obtained the same way as the model.The fitness of the solutionX is expressed as follows:(20)FSIX=13S1X−S12+S2X−S22+S3X−S32.


### 4.3.2. Experiment Results

The experiment is repeated 10 times for each algorithm. Table5 shows the best fitness, mean fitness, and standard deviation of test results on seismic inversion problem. The convergence curves and boxplots are shown in Figure 4. From the test results, the performance of OCDWOA for seismic inversion problem is the best.Table 5
Results of different algorithms.
ProblemAlgorithmBestMeanStdSeismic inversionWOA0.0071787540.01241271160.0038698430OWOA0.0061682620.01195090700.0053173170OEWOA0.0069015790.01120280110.0035417924IWOA0.0061365810.01289760650.0041582246OCDWOA0.0026622290.00559751190.0019518514The bold values mean the best results.Figure 4
Results of OCDWOA with respect to other algorithms. (a) Convergence graph of the algorithms. (b) Boxplots of the algorithms.
(a)(b)In Table5, the best fitness, mean fitness, and standard deviation of OCDWOA are all obviously smaller than those of the comparative algorithms. In Figure 4(a), the convergence rate of OCDWOA is better than those of the comparative algorithms. The convergence curve of OCDWOA is stable. In Figure 4(b), the stability of OCDWOA is the best among the comparative algorithms.


## 4.3.1. Seismic Inversion Problem

Seismic inversion is the key technology of quantitative reservoir prediction and prestack waveform inversion [44] is a newly developed seismic inversion technique. In prestack waveform inversion technique, the complex wave propagation effect is considered and the complete waveform information is fully used. The prestack waveform inversion method assumes that each CMP/CDP point of underground medium layer structure with local one-dimensional or level uses the reflectivity method forward prestack data [45]. The prestack waveform inversion method searches the best fitting with the actual prestack seismic data of underground through genetic algorithm, simulated annealing, and other global optimization algorithms [46]. At present, there are many problems in prestack waveform inversion technology that need to be further studied and the resolution of prestack waveform inversion technology needs to be improved. The existing optimization algorithms used in prestack waveform inversion are mainly traditional intelligent optimization algorithms. These algorithms are prone to premature convergence and local optimal solution when dealing with large data volume, highly nonlinear, multiparameter, and multiextremum prestack waveform inversion problems.In this paper, a theoretical geological model with 11 layers is established [47]. Each layer in the model contains 3 parameters: velocity of longitudinal wave (VP), velocity of shear wave (VS), and density (ρ). The values of model are shown in Table 4.Table 4
Values of the theoretical geological model.
VP (m·s−1)VS (m·s−1)ρ (g·cm−3)Layer 119808512Layer 2218516052.2Layer 319838562Layer 4220016222.2Layer 5243715102.3Layer 6266313802.4Layer 719838512Layer 8243715012.3Layer 919808562Layer 10265613802.4Layer 1119838512In this paper, the reflection coefficients of the model will be calculated by Zoeppritz equation as follows:(17)RPRSTPTS=−sinθ1−cosϕ1sinθ2cosϕ2cosθ1−sinϕ1cosθ2−sinϕ2sin2θ1VP1VS1cos2ϕ1ρ2VS22VP1ρ1VS12VP2cos2ϕ1ρ2VS2VP1ρ1VS12cos2ϕ2−cos2ϕ1VS1VP1sin2ϕ1ρ2VP2ρ1VP1cos2ϕ2−ρ2VS2ρ1VP1sin2ϕ2−1sinθ1cosθ1sin2θ1cos2ϕ1,where VP1, VS1, and ρ1 are the parameters of the up-layer, VP2, VS2, and ρ2 are the parameters of the low-layer, θ1 and θ2 are the incident angle and the refraction angle of the P-wave, respectively, and Φ1 and Φ2 are the incident angle and the refraction angle of the S-wave, respectively.Φ1, θ2, and Φ2 are calculated as follows, respectively:(18)ϕ1=arcsinθ1×VS1VP1,θ2=arcsinθ1×VS2VP1,ϕ2=arcsinθ1×VP2VP1.Incident anglesθ1 of 50 reflection coefficients are from 1° to 50°. The 50 seismic records are obtained by convolving the reflection coefficients with the Ricker wavelet with 35 Hz as follows:(19)S=R∗Rw,where S is the seismic record, R is the reflection coefficient, and Rw is the ricker wavelet (35 Hz).The 50 seismic records will be divided into three groups. The three groups are small incident angle group (Group 1), medium incident angle group (Group 2), and large incident angle group (Group 3). Group 1 consists of the seismic records whose incident angles are from 1° to 15°. Group 2 consists of the seismic records whose incident angles are from 16° to 32°. Group 3 consists of the seismic records whose incident angles are from 33° to 50°. Superimposed seismic records S1, S2, and S3 are obtained by synthesizing three groups of seismic records, respectively. The superimposed seismic records S1(X), S2(X), and S3(X) of the solution X are obtained the same way as the model.The fitness of the solutionX is expressed as follows:(20)FSIX=13S1X−S12+S2X−S22+S3X−S32.


## 4.3.2. Experiment Results

The experiment is repeated 10 times for each algorithm. Table5 shows the best fitness, mean fitness, and standard deviation of test results on seismic inversion problem. The convergence curves and boxplots are shown in Figure 4. From the test results, the performance of OCDWOA for seismic inversion problem is the best.Table 5
Results of different algorithms.
ProblemAlgorithmBestMeanStdSeismic inversionWOA0.0071787540.01241271160.0038698430OWOA0.0061682620.01195090700.0053173170OEWOA0.0069015790.01120280110.0035417924IWOA0.0061365810.01289760650.0041582246OCDWOA0.0026622290.00559751190.0019518514The bold values mean the best results.Figure 4
Results of OCDWOA with respect to other algorithms. (a) Convergence graph of the algorithms. (b) Boxplots of the algorithms.
(a)(b)In Table5, the best fitness, mean fitness, and standard deviation of OCDWOA are all obviously smaller than those of the comparative algorithms. In Figure 4(a), the convergence rate of OCDWOA is better than those of the comparative algorithms. The convergence curve of OCDWOA is stable. In Figure 4(b), the stability of OCDWOA is the best among the comparative algorithms.


## 5. Conclusions

In this paper, a modified WOA named OCDWOA is proposed. In order to balance the exploitation phase and exploration phase, a nonlinear parameter design is used. OBL, DE operator, and a density-based clustering strategy are introduced into the proposed algorithm to improve the performance of the search phase and the ability of jumping out of the local optimal solution. The proposed algorithm is tested on 19 optimization benchmark functions and a seismic inversion problem. The proposed algorithm is compared with WOA and the other three variants of WOA. Experimental results show that the proposed algorithm can enhance the performance of WOA and is better than other compared algorithms.Our future work will pay attention to the following two points: first, we would like to research the relation between the parameters and the performance of the OCDWOA. Second, we would like to apply the OCDWOA to machine learning.



---

*Source: 9159130-2022-01-24.xml*