Dynamic Harris Hawks Optimization with Mutation Mechanism for Satellite Image Segmentation

: In this paper, a novel satellite image segmentation technique based on dynamic Harris hawks optimization with a mutation mechanism (DHHO / M) is proposed. Compared with the original Harris hawks optimization (HHO), the dynamic control parameter strategy and mutation operator used in DHHO / M can avoid falling into the local optimum and e ﬃ ciently enhance the search capability. To evaluate the performance of the proposed method, a series of experiments are carried out on various satellite images. Eight advanced thresholding approaches are selected for comparison. Three criteria are adopted to determine the segmentation thresholds, namely Kapur’s entropy, Tsallis entropy, and Otsu between-class variance. Furthermore, four oil pollution images are used to further assess the practicality and feasibility of the proposed method on real engineering problem. The experimental results illustrate that the DHHO / M based thresholding technique is superior to others in the following three aspects: ﬁtness function evaluation, image segmentation e ﬀ ect, and statistical tests.


Introduction
Image segmentation is a fundamental and crucial stage in some applications, such as computer vision, pattern recognition, image classification, etc. [1][2][3][4].Generally speaking, the purpose of the segmentation operation is to partition the grayscale or color image into several non-overlapping regions with unique characteristics.The gray levels of the pixels in the same area are roughly the same, while the pixels from different areas are significantly different.In the last few decades, many segmentation methods have been proposed by scholars, such as clustering, fractal-wavelet modeling, region growing, and thresholding [5][6][7].Among all available techniques, threshold-based (thresholding) method is extensively used due its simplicity and efficiency [8,9].
Thresholding technique can be mainly divided into two categories: bi-level and multilevel [10].Bi-level thresholding splits the image into two classes according to the intensity.However, if the image contains much information and multiple objects, the bi-level thresholding is not competent.Therefore, as an extension of bi-level thresholding, multilevel thresholding is proposed to handle the current issues [11].For example, Díaz-Cortés et al. segment the breast thermograms using multilevel thresholding, which provides a highly reliable clinical decision support and advances the progress of thresholding approach in medical imaging domain [12].Also, multilevel thresholding is adopted by Bhandari et al. to process the satellite images, which improves the segmentation accuracy and obtains high-quality segmented images [13].In fact, researchers have introduced numerous approaches based on different criteria to select the segmentation thresholds over the past few decades, such as Kapur's entropy [14], Tsallis entropy [15], and Otsu between-class variance [16].Kapur's entropy measures the homogeneity of segmented regions by maximizing the histogram entropy.Tsallis entropy also determines the thresholds by maximizing the entropy, but there are certain constraints.Otsu's technique maximizes the between-class variance of the segmented image.
Multilevel thresholding can be considered as a NP-hard problem with high complexity [10].To be more specific, the computational complexity increases exponentially with the number of thresholds, which greatly affects the efficiency and feasibility.To overcome this drawback, meta-heuristic algorithms are introduced into this domain and combined with the thresholding technique.In 2017, He and Huang proposed a modified firefly algorithm (MFA) based method to select the segmentation thresholds [11].Kapur's entropy, Otsu between-class variance, and minimum cross entropy are served as fitness functions.Experimental results reveal that the proposed technique can efficiently reduce the computational complexity and produce segmented images with more details.Khairuzzaman and Chaudhury combined the grey wolf optimizer (GWO) with thresholding for grayscale image segmentation [7].Kapur's entropy and Otsu's method are used to select the segmentation thresholds.It can be found from the results that the proposed approach can accurately determine the segmentation thresholds and significantly reduce the computational complexity.In 2017, Oliva et al. presented an efficient magnetic resonance (MR) image segmentation method based on crow search algorithm (CSA) and minimum cross entropy [17].The results of the comparison experiment show that the well-delimited regions obtained are easier to distinguish than other techniques [17].Furthermore, some other algorithms or their modified versions have also been introduced into this domain, such as whale optimization algorithm (WOA) [18], multi-verse optimizer (MVO) [19], grasshopper optimization algorithm (GOA) [20], social spiders optimization (SSO) [21], krill herd optimization (KHO) [22], and cuckoo search (CS) [23] as well as Lévy flight firefly algorithm (LFA) [24], hybrid differential evolution (hjDE) [25], adaptive wind driven optimization (AWDO) [26], etc.These promising results motivate us to apply some other efficient meta-heuristic algorithms to multilevel image thresholding.Besides, as stated in No-Free Lunch Theorem (NFL) [27], there is no algorithm can solve all optimization problems.Thus, the recently proposed and uninvestigated algorithm also has potential.
Harris hawks optimization (HHO) is a high-performance, population-based, gradient-free optimization technique, which was proposed by Heidari et al. in 2019 [28].Inspired from the cooperative behavior and chasing style of Harris hawks in nature, the mathematical model of HHO was established [28].Similar to other meta-heuristic algorithms, the exploration and exploitation stages are also included in HHO.Each Harris hawk represents a candidate solution in the search space and the best solution obtained so far is considered as the prey or approximate optimal solution.During the whole process of predation, the Harris hawks patiently wait for the prey, besiege the prey to make it exhausted, and finally attack the prey with a surprise pounce.According to [28], HHO has presented better or occasionally competitive results than well-established techniques on 29 benchmark problems and several real-world engineering tasks.These phenomena illustrate the potential and the superiority of HHO, which also motivates us to apply this powerful algorithm to multilevel image thresholding.
No algorithm is perfect.Each algorithm requires some adjustments to fit the current problem, as does the HHO algorithm.For example, the exploration and exploitation phases of HHO are unbalanced.To be more specific, the absolute value of the parameter E cannot be larger than 1 during the latter half of the iteration (e.g., the last 250 of 500 maximum iterations), that is, the search agent does not perform a global search, although the current region may not be promising.For some complex optimization tasks, there is no guarantee that the population has gathered near the global optimum at the end of the exploration phase, resulting in premature convergence and local optimization.Therefore, two improvement strategies are introduced in this paper to enhance the optimization ability of the original HHO.As for the former, a novel dynamic control parameter strategy is presented to avoid trapping into the local optimum.The essence of the improvement strategy is to add a disturbance term to the update formulation of the escaping energy (E).Sine and cosine functions are used in combination to control where the disturbance peak appears.A Gaussian distribution is also adopted to increase randomness.In the experiment, not only the influence of parameter change on the performance of disturbance term, but also the reason for using the Gaussian distribution instead of the Levy distribution is analyzed.As for the latter, a famous mutation operator, known as DE/best/2 is utilized to improve the global search efficiency and population diversity [29].Then, the modified version of HHO, known as dynamic Harris hawks optimization with a mutation mechanism (DHHO/M) is applied to multilevel image thresholding.
Satellite imagery contains a wealth of information and is crucial to environmental resource monitoring and evaluation, such as forest cover, wetland resources, climate change, etc.Although remote sensing technology has many advantages (e.g., fast update cycle, fewer interference factors, and saving manpower and material resources), it is still a challenge to extract boundaries, locate objects and separate regions in high-resolution satellite images [30].To solve the above problems, DHHO/M is combined with thresholding technique for satellite image segmentation.A series of experiment are carried out to evaluate the performance of the proposed method.The experiments can be mainly divided into four parts.For the first part, the impact of both the improvement strategies on performance is investigated and analyzed.For the second and the most important part, the DHHO/M based approach is compared with four advanced thresholding techniques on satellite images.Various performance metrics have been taken into consideration, such as objective function value, standard deviation (Std), peak signal to noise ratio (PSNR), mean squared error (MSE), structural similarity index (SSIM), feature similarity index (FSIM), and convergence property as well as the Wilcoxon's rank sum test and Friedman test.For the third part, DHHO/M is compared with other thresholding techniques based on different criteria.For the fourth part, four oil pollution images are used to further assess its practicality and feasibility on real engineering problems.
The contributions of this paper have three aspects: 1.
Introduce the recently proposed HHO into multilevel image thresholding.To the best of the authors' knowledge, this attempt has not been made yet.

2.
Dynamic control parameter strategy and mutation mechanism are used to improve the search efficiency of the original HHO.

3.
Objectively and comprehensively evaluate the performance of the proposed technique.
The structure of this paper is described as follows: in Section 2, the problem statement and definitions of Kapur's entropy, Tsallis entropy, and Otsu's method are given.Section 3 briefly reviews the original HHO.In Section 4, the proposed DHHO/M based image thresholding technique is introduced in details.Section 3 presents the experimental results and relevant discussions.Finally, the conclusion and future research direction are represented in the last section.

Problem Statement
For multilevel image thresholding, the ultimate goal is to select the combination of segmentation thresholds.Assume that the given image is segmented into n + 1 classes {C 0 , C 1 , . . ., C n } by n thresholds {T 1 , T 2 , . . ., T n }, which can be defined as follows [31]: where I(i, j) denotes the gray level of the (i, j)th pixel, and L is the number of gray levels in the given image.
Furthermore, the combination of the optimal thresholds is determined by optimizing (maximizing or minimizing) the objective function, which is mathematically expressed as the following equations: where f represents the objective function and T * 1 , T * 2 , . . ., T * n is the combination of the optimal thresholds.

Kapur's Entropy
Kapur's entropy is a famous thresholding method that selects the thresholds based on the entropy of the segmented classes [14].The objective function can be defined as: where ω i = T i+1 −1 j=T i p j , p j = h( j)/N.ω i indicates the sum of the probabilities of the pixels in the ith class, p j is the proportion of pixels in each gray level to the total, h( j) represents the frequency, and N is the number of all pixels.Besides, for multilevel image thresholding, T 0 = 0 and T n+1 = L.

Tsallis Entropy
Tsallis entropy is extensively used by researchers in the field of image thresholding to choose the optimal combination of thresholds [15].The objective function is mathematically represented as: where q is a constant equal to 4, the value selected is the same as [30,32].

Otsu Between-Class Variance
Otsu's method was proposed in 1979, which maximizes the between-class variance of segmented classes to determine the thresholds [16].The objective function is presented as follows: where , µ = L−1 j=0 jp j .µ is the mean intensity for whole image.The mathematical form of above three optimization problems can be expressed as: Subject to x i must be integer, i = As discussed above, the computational complexity of the thresholding techniques will increase exponentially with the number of thresholds.Thus, meta-heuristic algorithms are used to optimize the objective functions in Equations ( 4)-( 6) to improve the efficiency.

Harris Hawks Optimization
HHO is a novel population-based, gradient-free optimization technique, which was proposed by Heidari et al. in 2019 [28].HHO simulates the predation, surprise pounce, and attacking behaviors of Harris hawks in nature.Similar to other meta-heuristic algorithms, HHO also includes two optimization stages, namely exploration and exploitation (see Figure 1), which are described in the following subsections.
As discussed above, the computational complexity of the thresholding techniques will increase exponentially with the number of thresholds.Thus, meta-heuristic algorithms are used to optimize the objective functions in Equations ( 4)-( 6) to improve the efficiency.

Harris Hawks Optimization
HHO is a novel population-based, gradient-free optimization technique, which was proposed by Heidari et al. in 2019 [28].HHO simulates the predation, surprise pounce, and attacking behaviors of Harris hawks in nature.Similar to other meta-heuristic algorithms, HHO also includes two optimization stages, namely exploration and exploitation (see Figure 1), which are described in the following subsections.

Exploration Stage
Harris hawks have insightful eyes that help them track and detect prey, however, sometimes the prey is not easy to find.Then, the Harris hawks will perch, wait patiently, and may last for hours.In HHO, above behaviors are modeled as the exploration stage, which can be expressed as: where  is the position of th individual in ( + 1)th iteration,  denotes the position of the rabbit (prey), and  is a random number in the interval [0,1] that converts the two strategies. ,  ,  , and  are also random numbers inside [0,1]. and  are the upper and lower bounds of the given optimization problem respectively. represents the average position of the population that can be calculated as follows: where  is the size of the population, and  is the position of th individual in th iteration.

Exploration Stage
Harris hawks have insightful eyes that help them track and detect prey, however, sometimes the prey is not easy to find.Then, the Harris hawks will perch, wait patiently, and may last for hours.In HHO, above behaviors are modeled as the exploration stage, which can be expressed as: where X t+1 i is the position of ith individual in (t + 1)th iteration, X rabbit denotes the position of the rabbit (prey), and q is a random number in the interval [0,1] that converts the two strategies.r 1 , r 2 , r 3 , and r 4 are also random numbers inside [0,1].LB and UB are the upper and lower bounds of the given optimization problem respectively.X t m represents the average position of the population that can be calculated as follows: where N is the size of the population, and X t i is the position of ith individual in tth iteration.

Transition from Exploration to Exploitation
The transition between exploration and exploitation is critical to the performance of meta-heuristic algorithms.In HHO, the escaping energy of the rabbit, known as E, is used to convert these two phases.The value of E decreases as the number of iterations increases, which can be mathematically modelled as: where E 0 is a random number that changes in the interval [−1,1], t denotes the current iteration, and t max represents the maximum number of iterations.More specifically, the value of E decreases from [−2,2] to 0 as the iterative process.If |E| ≥ 1, the exploration stage is used to search for the prey; If |E| < 1, the exploitation stage is adopted to exploit the promising area.

Exploitation Stage
After detecting the prey, the Harris hawks tend to attack it using surprise pounce.The actual predation process is often very complicated, the prey has a chance to escape and the Harris hawks also react differently according to the behavior of the prey.In order to better simulate the real situation, four strategies are used in the exploitation stage.A random number (r) is used to determine whether the prey has successfully escaped.The case r < 0.5 indicates successful escape, while r ≥ 0.5 shows unsuccessful escape.The escaping energy of the prey (E) affects the behavior of the Harris hawks.If |E| ≥ 0.5, the soft besiege occurs; if |E| < 0.5, the hard besiege happens [28].

A. Soft besiege
Considering r ≥ 0.5 and |E| ≥ 0.5, the rabbit has enough energy and keeps trying to escape, even though these acts are in vain and will only make the prey itself exhausted.At the same time, the Harris hawks softly encircle the rabbit and then attack it using surprise pounce.The mathematical model is given as follows: where ∆X t i indicates the difference between the position vector of the rabbit and the current individual.r 5 is a random number inside [0,1], and J = 2(1 − r 5 ) indicates the random jump strength of the rabbit during the escape process.

B. Hard besiege
Considering r ≥ 0.5 and |E| < 0.5, the rabbit becomes so exhausted and does not have enough escaping energy.Meanwhile, the Harris hawks hardly encircle the rabbit to conduct the final pounce.This process can be modeled as: C. Soft besiege with progressive rapid dives Considering r < 0.5 and |E| ≥ 0.5, the rabbit has enough energy and managed to escape.In the meantime, a more intelligent soft besiege is performed before the surprise pounce.Note that the position update of the Harris hawks in this situation is a two-step process.If the first step is not an improved move, then the second step will be executed.For the former part, it can be mathematically modelled as: For the latter part, the Levy flight is used to simulate the rapid, abrupt, and irregular movement of Harris hawks when chasing the rabbit.The formulation for position update is presented as: where dim is the dimensions of the optimization problem, S is a random vector of size 1 × dim, and Levy denotes the Levy distribution that is defined as: where u and v are random values in the interval [0,1].β is a constant equal to 1.5.
Then, for a problem to be minimized, the mathematical model of the whole process at this stage is: where f is the fitness function for the given optimization problem.

D. Hard besiege with progressive rapid dives
Considering r < 0.5 and |E| < 0.5, the rabbit does not have enough energy but it successfully escapes.In this situation, the Harris hawks perform a hard besiege before the surprise pounce, and they try to bring the whole group closer to the prey, not just the individual.Therefore, the formulation for position update of the Harris hawks is given as follows: where Furthermore, the pseudocode of HHO is represented in Algorithm 1.

Algorithm 1 Pseudocode of Harris hawks optimization (HHO) algorithm
Input: The size of population N, maximum number of iterations t max .
Output: The position of the rabbit and the corresponding fitness function value.

1.
Initialize the position of the hawks X i and the rabbit X rabbit .

2.
Initialize the fitness values of the hawks f i and the rabbit f rabbit .

3.
Set the dimensions of the optimization problem dim.
Check the boundary and evaluate the fitness value of each hawk f i .6.
Update the location X rabbit and fitness value f rabbit of rabbit if there is a better one.7.
Update the energy of the rabbit E using Equation (9).9.

Proposed Dynamic Harris Hawks Optimization with Mutation Mechanism
In this section, the proposed DHHO/M algorithm based image segmentation technique is introduced in detail.For multilevel color image segmentation problem, the input to the method was the histogram of each color component, and the output was the segmentation thresholds.This process was implemented by optimizing the objective functions in Equations ( 4)-( 6) using the DHHO/M algorithm.The following subsections will discuss the two improvement strategies respectively, and then give the algorithm steps.Finally, the time complexity of DHHO/M and HHO is analyzed.

Dynamic Control Parameter Strategy
In the original HHO algorithm, the transition of the exploration and exploitation phase is based on the escaping energy (E) of the rabbit.More specifically, the absolute value of the parameter E decreased from 2 to 0 throughout the whole iteration.If |E| > 1, HHO explored the search space to determine the promising area.On the contrary, if |E| < 1, the exploitation strategy was adopted to enhance the local search efficiency.However, the main drawback of this transition process is that the absolute value of the parameter E cannot be larger than 1 during the latter half of the iteration (e.g., the last 250 of 500 maximum iterations), that is, the search agents do not perform global search during the second half of the iteration, although the current region may not be promising.For some complex optimization tasks, there was no guarantee that the population has gathered near the global optimum at the end of the exploration phase, resulting in premature convergence and local optimal values.Therefore, a novel dynamic control parameter strategy is introduced with the purpose of jumping out of the local optimum under the premise of ensuring accuracy.The essence of the improvement strategy is to add a disturbance term to the update formulation of E in Equation ( 9).The disturbance term and modified update formulation of E are presented as follows: where sin and cos denote the sine and cosine functions respectively.randn is a random number subject to Gaussian distribution.t and t max represent the current iteration and the maximum number of iterations respectively.α is a constant that determines where the disturbance peak appears.

292
In order to more intuitively illustrate the impact of dynamic control parameter strategy on the 293 escaping energy of the rabbit, the variation of the parameter  with the number of iterations is also 294 given in Figure 3.It is worth mentioning that the values of  was not only the key to global and local 295 phases transition, but also affects the convergence of the algorithm because it was part of the position 296 update formulation.To be more specific, a larger value indicates more exploration, while a smaller 297 value indicates more exploitation.As can be observed from the figure, when  = 2, the absolute 298 value of  in the late iteration can be abruptly greater than 1, but the amplitude of the disturbance 299 term is too large, which seriously affects the convergence performance.In the case of  = 3, the 300 In order to more intuitively illustrate the impact of dynamic control parameter strategy on the escaping energy of the rabbit, the variation of the parameter E with the number of iterations is also given in Figure 3.It is worth mentioning that the values of E was not only the key to global and local phases transition, but also affects the convergence of the algorithm because it was part of the position update formulation.To be more specific, a larger value indicates more exploration, while a smaller value indicates more exploitation.As can be observed from the figure, when α = 2, the absolute value of E in the late iteration can be abruptly greater than 1, but the amplitude of the disturbance term is too large, which seriously affects the convergence performance.In the case of α = 3, the amplitude of the perturbation term is too small, which cannot enhance the ability of jumping out of the local optimum, but will not affect the performance of the algorithm itself.Considering α = 2.5, it can be found that the convergence property was not greatly affected, and the global search strategy also had an opportunity to be adopted in the later iteration.As far as the authors are concerned, for optimizing some fitness functions that contain local best, setting α to 3 may reduce the accuracy; setting α to 2 made the results obtained similar to the original algorithm; while setting α to 2.5 can improve the performance of the algorithm.These hypotheses will be verified in the experimental section.In order to more intuitively illustrate the impact of dynamic control parameter strategy on the escaping energy of the rabbit, the variation of the parameter  with the number of iterations is also given in Figure 3.It is worth mentioning that the values of  was not only the key to global and local phases transition, but also affects the convergence of the algorithm because it was part of the position update formulation.To be more specific, a larger value indicates more exploration, while a smaller value indicates more exploitation.As can be observed from the figure, when  = 2, the absolute value of  in the late iteration can be abruptly greater than 1, but the amplitude of the disturbance term is too large, which seriously affects the convergence performance.In the case of  = 3, the amplitude of the perturbation term is too small, which cannot enhance the ability of jumping out of the local optimum, but will not affect the performance of the algorithm itself.Considering  = 2.5, it can be found that the convergence property was not greatly affected, and the global search strategy also had an opportunity to be adopted in the later iteration.As far as the authors are concerned, for optimizing some fitness functions that contain local best, setting  to 3 may reduce the accuracy; setting  to 2 made the results obtained similar to the original algorithm; while setting  to 2.5 can improve the performance of the algorithm.These hypotheses will be verified in the experimental section.In addition, it was necessary to explain the reason why the Gaussian distribution is used instead of the Levy distribution in the disturbance term.Although the characteristic of usually small steps and occasional large jumps in Levy distribution fits well with the idea of the improvement strategy, sometimes there are some extreme cases as shown in Figure 4, which can seriously affect the performance of the algorithm.More specifically, the disturbance of uncertain position makes the population evolve slowly and even degenerate.
Reducing the step size of Levy flight can theoretically improve this situation, but in fact, it is difficult to determine an appropriate value to meet the desired requirements.In our study, if the step size is larger than 0.05, disturbances with large amplitude will always appear; if the step size is smaller than 0.01, the disturbance is not enough to affect the value of E; while the value between the two is still not ideal.Therefore, the Gaussian distribution is more suitable for the current problem, and the amplitude generated by it is acceptable compared with the Levy distribution.In addition, it was necessary to explain the reason why the Gaussian distribution is used instead 312 of the Levy distribution in the disturbance term.Although the characteristic of usually small steps 313 and occasional large jumps in Levy distribution fits well with the idea of the improvement strategy, 314 sometimes there are some extreme cases as shown in Figure 4, which can seriously affect the 315 performance of the algorithm.More specifically, the disturbance of uncertain position makes the 316 population evolve slowly and even degenerate.317 Reducing the step size of Levy flight can theoretically improve this situation, but in fact, it is 318 difficult to determine an appropriate value to meet the desired requirements.In our study, if the step 319 size is larger than 0.05, disturbances with large amplitude will always appear; if the step size is 320 smaller than 0.01, the disturbance is not enough to affect the value of ; while the value between the 321 two is still not ideal.Therefore, the Gaussian distribution is more suitable for the current problem, 322 and the amplitude generated by it is acceptable compared with the Levy distribution.323 324  Furthermore, according to our experimental statistics, the disturbance peaks of our proposed strategy usually appearred during the 300th to 400th iterations (the maximum number of iterations was 500), especially the 350th to the 400th iterations.

Mutation Mechanism
In order to enhance the global search efficiency of HHO, mutation mechanism was introduced into the exploration stage.The famous DE/best/2 mutation operator was adopted to replace the original location update strategy in Equation ( 7) [29].The formulation is given as follows: where x t+1 i indicates the position of the ith individual at the tth iteration, and F is the scaling factor.r 1 , r 2 , r 3 , and r 4 are different integers selected from the range [1, N] and distinct from i. N denotes the size of population.Thus, Equation ( 7) can be rewritten as follows: where rand1 and rand2 are random numbers in the interval [0,1].The operators of differential evolution (DE) have been widely used as part of modified or hybrid algorithm to improve the optimization capability and search efficiency.This is also the motivation for us to adopt the mutation operator mentioned above.For example, Xiong et al. proposed a hybrid algorithm of whale optimization algorithm (WOA) and DE, known as DE/WOA, for extracting parameters of solar photovoltaic models [33].DE/best/2 is utilized in this proposed algorithm with the purpose of improving the ability of exploring the search space and locating the region of global optimum [33].Also, Jadon et al. used DE to modify the onlooker bee phase of artificial bee colony (ABC) algorithm [34].DE/best/1 is adopted to accelerate the convergence speed and maintain population diversity.Furthermore, DE has also been embedded or merged into some other improved versions of meta-heuristic algorithm, such as moth search algorithm based DE (MSDE) [35], improved whale optimization algorithm (IWOA) [36], self-adaptive gravitational search algorithm and DE (SGSADE) [37], improved dragonfly algorithm (IDA) [38], gbest-guided ABC (GABC) [39], hybrid whale optimization algorithm with DE (WOA-DE) [18], chaotic opposition-based grey wolf optimization algorithm based on DE and disruption operator (COGWO2D) [40], etc.
In fact, there are other types of mutation and crossover operators that have been proposed by researchers to enhance the optimization ability.In 2019, Xu et al. introduced a series of new variants of moth-flame optimization (MFO) by combining MFO with Gaussian mutation (GM), Cauchy mutation (CM), Lévy mutation (LM) or the combination of GM, CM and LM [41].The experimental results demonstrated that the three strategies can significantly boost exploration and exploitation capabilities of the basic MFO.For the crossover operators, Kita et al. analyzed the unimodal normal distribution crossover (UNDX) by discussing the importance of the distribution and statistics of the offspring yielded by a crossover operator for its evaluation [42].Akimoto et al. proposed an adaptive real-coded ensemble crossover (AREX) that combines the adaptation of expansion rate technique and the crossover mean descent technique in 2009 [43].Experimental results showed that the proposed algorithm enlarges the classes of functions that original methods can solve.In 2015, Ariyarit and Kanazaki introduced a modified version of genetic algorithm (GA) with multi-modal distribution crossover (MMDX) [44].Blend crossover (BLX) [45] and UNDX were used as comparison algorithms.It can be found from the experimental results that MMDX is more competitive and can maintain higher population diversity.
In addition, experiments are performed on eight satellite images (see Figure 5) to determine the optimal parameter settings of the proposed algorithm.The sum of ranks with different combination of α and F is given in Table 1.The lower the ranking value, the better the performance.As can be observed, α = 2.5 and F = 0.5 present the best results in most cases, thus this combination of parameter values is adopted in the following experiments.28.Set population size N and maximum number of iterations t max .29. Set the dimensions of the optimization problem dim, namely the number of thresholds.30.While (termination condition is not met (t < t max )) do 31.Check the boundary and evaluate the fitness value of each hawk f i using Equations ( 4)- (6).32.Update the location X rabbit and fitness value f rabbit of rabbit if there is a better one.33.
Update the energy of the rabbit E using Equation (21).% Dynamic control parameter strategy 35.

48.
End For 49. End While Fitness function (Kapur's entropy) Input: Histogram of a color component, and segmentation thresholds X i .Output: Fitness function value f i .

1.
The histogram is divided into n + 1 parts by n thresholds.

2.
Calculate the proportion of pixels in each gray level p j , j ∈ [0, 255] to the total based on the histogram.

3.
Compute the sum of the probabilities of the pixels (ω k , k ∈ [0, n]) contained in each part.

5.
The sum of the entropies of all parts represents the fitness function value.6.

Computational Complexity
The computational complexity of DHHO/M and HHO algorithms mainly depends on three processes: initialization, function evaluation, and position update.The improvement strategy of this paper only replaces the relevant formulations, and does not add additional judgments, loops or other commands.Therefore, the computational complexity of DHHO/M is approximately the same as that of HHO, which can be computed as follows: where N is the size of population, Co f is the cost of function evaluation, dim is the dimension of the given optimization problem, and t max is the maximum number of the iterations [46].
In fact, the proposed DHHO/M runs slightly faster than the original HHO, which will be demonstrated and analyzed in subsequent experiments.

Experimental Setup and Database
In order to assess the performance of the proposed approach, a series of experiments were carried out on various images.More specifically, the experiments can be divided into four parts.For the first part of the experiment, the performance influence of each strategy (dynamic control parameter and mutation mechanism) is investigated, and the relevant algorithms are denoted by DHHO and HHO/M, respectively.Two DHHO variants with different parameter values were also included to further evaluate the effect of the parameter (α).The algorithm with parameter α equal to three is represented as DHHO + , and DHHO -indicates α = 2.For the second part of the experiment, the proposed method is compared with the traditional HHO and four advanced multilevel thresholding techniques, namely teaching-learning-based optimization (TLBO) [47], whale optimization algorithm with TH heuristic (WOA-TH) [48], improved differential search algorithm (IDSA) [49], and beta differential algorithm (BDE) [50].The experiment is conducted on eight satellite images (see Figure 5) with low and high threshold levels [51].For the third part of the experiment, DHHO/M is compared with other thresholding techniques based on different criteria to illustrate its feasibility and universality, such as Tsallis entropy based modified grasshopper optimization algorithm (MGOA) [20] and modified artificial bee colony (MABC) [13], as well as Otsu method based modified flower pollination algorithm (MFPA) [48] and grey wolf optimizer (GWO) [7].For the fourth part of the experiment, four oil pollution images are used to further evaluate its practicality and feasibility on real engineering problems.
It is worth noting that the parameters of each thresholding method used are the same as the original literature, except for the population size N set to 30 and the number of iterations t max set to 500 for fair comparison.The parameter setting can be found in Table 2. Besides, the initial population of all algorithms is randomly generated in the solution space, and the termination condition in this work is to reach the maximum number of iterations t max .A number of 50 independent runs are conducted, and the best result is highlighted in boldface.All experiments are executed on a computer with "Microsoft Windows 10" system and "8GB" memory space under Matlab2017.

Performance Metrics.
Several performance metrics are briefed in this section, which can be observed in Table 3.Among the available measures, average fitness function value and standard deviation (Std) are used to evaluate the optimization capability of meta-heuristic algorithm, while peak signal to noise ratio (PSNR), mean squared error (MSE), structural similarity index (SSIM), and feature similarity index (FSIM) are utilized to assess the quality of the segmented image.Furthermore, Wilcoxon's rank sum test and Friedman test were adopted for statistical analysis.

No. Measures Formulation Reference
Respective energy Respective energy

Performance Metrics.
Several performance metrics are briefed in this section, which can be observed in Table 3.Among the available measures, average fitness function value and standard deviation (Std) are used to evaluate the optimization capability of meta-heuristic algorithm, while peak signal to noise ratio (PSNR), mean squared error (MSE), structural similarity index (SSIM), and feature similarity index (FSIM) are utilized to assess the quality of the segmented image.Furthermore, Wilcoxon's rank sum test and Friedman test were adopted for statistical analysis. 9.
Friedman test

Experimental Series 1: Performance on Mathematical Functions
In this section, the proposed algorithm is evaluated on three mathematical test functions with different characteristics, such as unimodal, multi-modal, and composite.The test functions are taken from the CEC2005 special session, which are known as F1, F9, and F14.It can be found from Table 4 (bold is the best) that the proposed method outperforms others in most cases, showing its remarkable performance on mathematical test functions.This promising result motivates us to apply it to multilevel image thresholding domain.

Experimental Series 2: Influence of Dynamic Control Parameter Strategy and Mutation Mechanism
In this section, the influence of each improvement strategy (dynamic control parameter and mutation mechanism) is investigated, and the results can be found in Table 5 (bold is the best).As can be observed, DHHO/M presents the highest fitness function value in most cases, while HHO/M and DHHO are ranked second and third respectively.This promising result indicates that introduction of each strategy can improve the performance of the original algorithm, and the combination of them works better.Specifically speaking, the dynamic control parameter strategy prevents falling into local optimum, and the mutation mechanism effectively enhances the global search capability.Thus, DHHO/M exhibits superior performance compared to other algorithms.
Another comparison is made between DHHO, HHO, DHHO + , and DHHO -.From the table, these algorithms can be arranged as DHHO > HHO ~DHHO + > DHHO -, which strongly supports the discussion in the previous section.Smaller value of parameter (α = 2) makes the amplitude of the disturbance term too large, although the diversity of the population is increased but the algorithm convergence ability is seriously affected.On the contrary, larger value of parameter (α = 3) produces a small disturbance.Although the local optimum cannot be avoided, it does not affect the accuracy of the algorithm, so the result is similar to the original algorithm.However, the value of parameter α of DHHO is between the above two, which can enhance the ability to jump out of local optimum while ensuring convergence property, and obtain the best result in this comparison.
Furthermore, the Std values obtained by each algorithm are also given in Table 5.It can be found that DHHO/M presents the lowest value in general, which shows better stability than other algorithms.More detailed and comprehensive stability testing and analysis can be found in the next subsection.

Experimental Series 3: Comparison with Other Advanced Methods on Satellite Images
In this section, the proposed method is compared with other advanced multilevel image segmentation approaches on satellite images.The performance of each algorithm is comprehensively and objectively assessed in terms of precision, stability, statistical significance, scalability (high dimension), and convergence property.

Segmentation Accuracy
Accuracy is one of the main indicators for measuring image segmentation technique.In this experiment, the precision of each method is evaluated from five aspects, namely average fitness function value, MSE, PSNR, SSIM, and FSIM values.Firstly, the threshold values and segmented images obtained are given in Table 6 and Figure 6, respectively.The threshold levels selected (K = 2, 3, 4, and 5) are the same as [47,49].As can be observed from Table 6, all methods give similar results when the threshold level is small; while the threshold values determined become different when K = 4 and 5.The reason for this phenomenon is that the difficulty of the current optimization problem increases with the number of thresholds, and some meta-heuristic algorithms with poor search ability are not competent, such as "Image3" at three threshold levels, "Image5" at three threshold levels, and "Image6" at five threshold levels.Considering the segmented images presented in Figure 6, it can be found that the images with high levels contain more information and details than that with low levels, because the of entropy of a given image reflects its average information content [18].Secondly, the average fitness function values obtained by all algorithms are given in Table 7 (bold is the best), which reflects the optimization capability of each meta-heuristic algorithm.Higher Kapur's entropy value indicates stronger search ability and higher precision.It can be observed that the values of all algorithms increase with the number of threshold levels, and DHHO/M outperforms other compared methods in general.This results strongly proves that high-quality segmented image with more information can be obtained at higher threshold level, such as K = 5.In DHHO/M, the introduction of dynamic control parameter strategy and mutation mechanism can improve the search efficiency of the algorithm and avoid premature convergence.Therefore, the superiority of the proposed algorithm is holistic rather than on a single image.
Thirdly, several performance metrics are introduced to assess the similarity between the segmented image and the original image, such as PSNR, MSE, SSIM, and FSIM.The experimental results are given in two tables.In Table 8, it can be found that the PSNR value increases with the number of threshold values, while MSE value is the opposite (bold is the best).Similarly, SSIM and FSIM values also increase with the number of threshold levels, as can be seen from Table 9 (bold is the best).These promising phenomena indicate that the quality of segmented image is improving gradually; while the results obtained by the proposed method are better than or at least comparable to other approaches, showing better feature and structure preserving ability.

Statistical Test
The experiments performed by each algorithm are the same, so the statistical test is necessary.Parametric statistical test is based on various assumptions, such as independence, normality, and synchronization [56].However, these assumptions probably not valid when analyzing meta-heuristic algorithms with stochastic property.Therefore, two famous nonparametric statistical tests are adopted in this study, namely, Wilcoxon's rank sum test [54] and Friedman test [55].As a brief review, the null hypothesis (H 0 ) assumes that there is no difference between the two compared algorithms (Wilcoxon's rank sum test) or the whole (Friedman test); while the alternative hypothesis H 1 indicates the difference.More detailed discussion of these two statistical tests can be found in [55].
The experimental results are given in three tables.Table 10 presents the p-values and h-values of Wilcoxon's rank sum test obtained by each comparison.It can be found that DHHO/M gives better results in 29 out of 32 cases (eight images and four threshold levels) for HHO, 30 cases for TLBO, 31 cases for WOA-TH, 30 cases for IDSA, and 30 cases for BDE, which shows significant difference between the proposed method and other approaches.As for the average rank of Friedman test given in Table 11, the proposed method ranks first in all cases, and the value of rank becomes smaller as the number of threshold level increases, indicating greater advantages over other compared methods (bold is the best).Table 12 gives the Chi-square (χ 2 ) value and p-value of Friedman test at different threshold levels.According to the Chi-square distribution table, the critical value for 5 (6 algorithm -1) degrees of freedom at 0.05 significant level is 11.07 [57].It can be observed from the table that the Chi-square values obtained are much larger than the critical value and the p-values obtained are much smaller than the significant level.The results illustrate that the advantage of the proposed algorithm becomes more obvious as the dimension of the optimization problem increases, and there is a significant difference between the available methods.

Computational Time
Computational time was also one of the important indicators to measure the efficiency of the algorithm.The experimental results at different threshold levels can be found in Table 13 (bold is the best).Obviously, the IDSA algorithm is the fastest, although it cannot give the best results in most cases.DHHO/M is a little faster than the traditional HHO, because the greedy selection used in exploitation stage consumes a lot of time; while the dynamic control parameter strategy enables the exploration phase can be adopted in the later iteration, which reduces the complexity to some extent.It is worth noting that evaluating a meta-heuristic algorithm should be comprehensive and objective, not just from a single aspect.The proposed method presents better or at least competitive results in terms of average fitness value, PSNR, MSE, SSIM, and FSIM, despite not the fastest.Thus, DHHO/M can be considered as an efficient technique for multilevel color image segmentation.In this section, all available thresholding approaches are compared at K=10, 15, and 20 to assess the performance on high dimensional tasks.Threshold levels selected are the same as [57].The average fitness function values obtained are presented in Table 14 (bold is the best).It can be observed that DHHO/M gives the best results in 22 out of 24 cases (eight images and three thresholds); TLBO and WOA-TH only give the best results in one case respectively, while other algorithms give none.As discussed in above section, the proposed method has shown certain advantages when K = 5, but it is more obvious in higher dimensions.The reason for this phenomenon is that each image can be considered as a different optimization problem, and there is no algorithm can handle all tasks [27].Even the meta-heuristic algorithm with strong search capability cannot show remarkable performance on all images.Therefore, DHHO/M is more competitive in the field of image segmentation than other algorithms, because it gives the best results in most cases.

Stability Analysis
The stability of each algorithm is analyzed and discussed in this section, including two performance measures, namely Std value and boxplot.The Std values obtained on satellite images can be found in Table 7, respectively.A lower Std value indicates better performance.As can be observed, the proposed approach presents the lowest Std values in most cases, showing remarkable stability, continuity and consistency.Furthermore, the boxplot with 20 threshold levels are given in Figure 7, because the difference between the algorithms is more obvious in high dimensions.It can be found that DHHO/M outperforms other approaches again.More specifically, the boxplot obtained is higher in position and more compact, indicating better data consistency.Compared with other algorithms, the proposed method has no bad points (in our study), which illustrates the superiority of DHHO/M based thresholding technique in satellite image segmentation.In this section, the convergence property of DHHO/M algorithm is evaluated.The convergence curves obtained by all algorithms at 20 threshold levels are presented in Figure 8.It is worth mentioning that these figures are drawn in a semi-logarithmic coordinate system for ease of observation.As can be seen, WOA-TH and IDSA are prone to premature convergence and local optimization, which is reflected in the shape of the curve.For example, under the circumstance of "Image3" segmentation, the fitness function values obtained by these two algorithms no longer change after 100 iterations, and obviously they do not get the best value at the end of iteration.Although the fast convergence speed is not a bad characteristic of meta-heuristic algorithm, local optimization seriously affects the accuracy as well as the quality of segmented image.On the contrary, the propose method can better balance the exploration and exploitation stages, neither premature convergence nor slow convergence.Therefore, the DHHO/M based technique can exhibit higher search efficiency in all selected satellite images.

Convergence Property
In this section, the convergence property of DHHO/M algorithm is evaluated.The convergence curves obtained by all algorithms at 20 threshold levels are presented in Figure 8.It is worth mentioning that these figures are drawn in a semi-logarithmic coordinate system for ease of observation.As can be seen, WOA-TH and IDSA are prone to premature convergence and local optimization, which is reflected in the shape of the curve.For example, under the circumstance of "Image3" segmentation, the fitness function values obtained by these two algorithms no longer change after 100 iterations, and obviously they do not get the best value at the end of iteration.Although the fast convergence speed is not a bad characteristic of meta-heuristic algorithm, local optimization seriously affects the accuracy as well as the quality of segmented image.On the contrary, the propose method can better balance the exploration and exploitation stages, neither premature convergence nor slow convergence.Therefore, the DHHO/M based technique can exhibit higher search efficiency in all selected satellite images.

Experimental Series 4: Performance Using Different Objective Functions
In this section, Tsallis entropy and Otsu between-class variance are served as objective function to assess the feasibility of DHHO/M.Eight different color images are selected from Figure 5.The experiment are conducted on high threshold levels and can be divided into two parts.For the former part, Tsallis entropy based thresholding techniques are selected for comparison, namely the MGOA and MABC.For the latter part, Otsu between-class variance based thresholding approaches are chosen for testing, namely the MFPA and GWO.
The experimental results of Tsallis entropy and Otsu's method are presented in Table 15 and Table 16, respectively (bold is the best).As can be observed, DHHO/M gives the best objective function values in all cases, indicating the strong optimization capability.Considering other metrics, the proposed algorithm also produces better or at least competitive results, which shows that the segmented images with more information and details can be obtained.To sum up, the superior performance of the proposed algorithm is not limited to the objective function adopted in this paper, but has potential in other image fields.

Experimental Series 4: Performance Using Different Objective Functions
In this section, Tsallis entropy and Otsu between-class variance are served as objective function to assess the feasibility of DHHO/M.Eight different color images are selected from Figure 5.The experiment are conducted on high threshold levels and can be divided into two parts.For the former part, Tsallis entropy based thresholding techniques are selected for comparison, namely the MGOA and MABC.For the latter part, Otsu between-class variance based thresholding approaches are chosen for testing, namely the MFPA and GWO.
The experimental results of Tsallis entropy and Otsu's method are presented in Tables 15 and 16, respectively (bold is the best).As can be observed, DHHO/M gives the best objective function values in all cases, indicating the strong optimization capability.Considering other metrics, the proposed algorithm also produces better or at least competitive results, which shows that the segmented images with more information and details can be obtained.To sum up, the superior performance of the proposed algorithm is not limited to the objective function adopted in this paper, but has potential in other image fields.In order to further verify the effectiveness of the proposed algorithm in solving practical engineering problems, four oil pollution images are selected for the experiment.These images were taken by the drone in the area of the eighth oil production plant, which can be found in Figure 9 [58].As can be observed, the oil pollution area of (b) is relatively obvious; while the remaining three images all have strong interference, especially (d).This phenomenon undoubtedly increases the difficulty of segmentation operation, which can be considered as a challenging engineering problem.

Experimental Series 5: Oil Pollution Image Segmentation
In order to further verify the effectiveness of the proposed algorithm in solving practical engineering problems, four oil pollution images are selected for the experiment.These images were taken by the drone in the area of the eighth oil production plant, which can be found in Figure 9 [58].As can be observed, the oil pollution area of (b) is relatively obvious; while the remaining three images all have strong interference, especially (d).This phenomenon undoubtedly increases the difficulty of segmentation operation, which can be considered as a challenging engineering problem.The segmented results of the proposed method are presented in Figure 10.Because there is no absolute standard for a real engineering problem, the authors manually labeled the target region and separated it.Then took it as the ground truth for experimental comparison.It can be found from the figures that the oil pollution in (a)-(c) has been successfully separated from the complex background, which is similar to the ground truth.Considering figure (d), the ideal oil pollution area consists of two parts, while the DHHO/M based thresholding technique incorrectly classifies them as a whole.But in fact, it has achieved the desired goal, because the main area of oil pollution has been efficiently identified.To sum up, although the proposed method cannot achieve perfect results for complex images with mixed background and strong interference, it is competent for most cases and can still be considered as a competitive technique for the segmentation of oil pollution image.
figures that the oil pollution in (a)-(c) has been successfully separated from the complex background, which is similar to the ground truth.Considering figure (d), the ideal oil pollution area consists of two parts, while the DHHO/M based thresholding technique incorrectly classifies them as a whole.But in fact, it has achieved the desired goal, because the main area of oil pollution has been efficiently identified.To sum up, although the proposed method cannot achieve perfect results for complex images with mixed background and strong interference, it is competent for most cases and can still be considered as a competitive technique for the segmentation of oil pollution image.

Conclusions
An efficient satellite image segmentation technique based on DHHO/M is proposed in this paper.Dynamic control parameter strategy and mutation mechanism are used to avoid trapping into the local optimum and improve the search efficiency.In order to validate the superiority of the proposed method, a series of experiments are conducted on various images.For the first part of the experiment, the experimental results indicate that both the improvement strategies adopted can

Conclusions
An efficient satellite image segmentation technique based on DHHO/M is proposed in this paper.Dynamic control parameter strategy and mutation mechanism are used to avoid trapping into the local optimum and improve the search efficiency.In order to validate the superiority of the proposed method, a series of experiments are conducted on various images.For the first part of the experiment, the experimental results indicate that both the improvement strategies adopted can enhance the optimization capability.For the second part of the experiment (the most crucial part), it can be observed that the DHHO/M based technique gives better results for satellite images in terms of objective function value, Std, PSNR, SSIM, FSIM, and convergence property as well as the Wilcoxon's rank sum test and Friedman test.For the third part of the experiment, the robustness of the proposed approach is assessed by the segmentation based on other criteria.For the last part of the experiment, the proposed method is applied to a real engineering problem, namely the segmentation of oil pollution image to further evaluate its practicality and feasibility.The experimental results strongly illustrate the remarkable performance of the DHHO/M based thresholding technique.
In the future, the authors will introduce more efficient technique to handle the images with mixed background and strong interference.Furthermore, due to the conflicts between different criteria, multi-objective optimization for image thresholding is also the research direction.

Figure 2 287 Figure 2 Figure 2 .
Figure 2 illustrates the variation of the disturbance term  under different values of parameter 288 .It can be found that the amplitude decreases as the value of the parameter  increases, and the 289 high amplitude disturbance occurs earlier when  = 2. 290

Figure 2 .
Figure 2. Variation of the disturbance term δ under different values of parameter α.

Figure 2 .
Figure 2. Variation of the disturbance term  under different values of parameter .

3 Figure 3 .Figure 3 .
Figure 3. Variation of the escaping energy of the rabbit  under different values of parameter .

Figure 4 .
Figure 4. Variation of the escaping energy of the rabbit  using Levy distribution

Figure 4 .
Figure 4. Variation of the escaping energy of the rabbit E using Levy distribution.

Figure 5 .
Figure 5. Original satellite images and the corresponding histograms for each of color channels (red, green and blue).

Figure 5 .
Figure 5. Original satellite images and the corresponding histograms for each of color channels (red, green and blue).

Figure 6 .
Figure 6.The segmented results of satellite images obtained by different algorithms.

Figure 6 .
Figure 6.The segmented results of satellite images obtained by different algorithms.

Figure 7 .
Figure 7.The boxplot of each method on satellite images (K = 20).

Figure 7 .
Figure 7.The boxplot of each method on satellite images (K = 20).

Figure 8 .
Figure 8.The convergence curves for fitness function on satellite images (K = 20).

Figure 8 .
Figure 8.The convergence curves for fitness function on satellite images (K = 20).

Figure 9 .
Figure 9. Original oil pollution images and the corresponding histograms for each of color channels (red, green and blue).

Figure 9 .
Figure 9. Original oil pollution images and the corresponding histograms for each of color channels (red, green and blue).

Figure 10 .
Figure 10.The segmented results and ground truth of the oil pollution images.Panels (a-d) indicate the segmented images, (e-h) show the ground truth, (i-l) indicate the classification result.

Figure 10 .
Figure 10.The segmented results and ground truth of the oil pollution images.Panels (a-d) indicate the segmented images, (e-h) show the ground truth, (i-l) indicate the classification result.

Table 1 .
Sum of ranks obtained by the proposed method (K = 5) under different of parameter values.

Table 2 .
Parameters of the algorithms.conducted,and the best result is highlighted in boldface.All experiments are executed on a computer 416 with "Microsoft Windows 10" system and "8GB" memory space under Matlab2017.417

Table 2 .
Parameters of the algorithms.

Table 3 .
Performance measures of the multilevel image segmentation methods.

Table 3 .
Performance measures of the multilevel image segmentation methods.

Table 4 .
The experimental results of different algorithms on mathematical test functions.

Table 5 .
The experimental results of HHO variants on satellite images.

Table 6 .
Comparison of optimal thresholds for different algorithms using Kapur's entropy at 2-5 levels.

Table 7 .
The average fitness values and std values of different algorithms at 2, 3, 4, and 5 levels.

Table 8 .
The PSNR values and MSE values of different algorithms at 2, 3, 4, and 5 levels.

Table 10 .
Results of Wilcoxon's rank sum test over all available satellite images at 2, 3, 4, and 5 levels.

Table 11 .
Friedman rank of different algorithms considering all experimental data at 2, 3, 4, and 5 levels.

Table 12 .
Results of Friedman test over all available satellite images at 2, 3, 4, and 5 levels.

Table 13 .
The average computation time (in seconds) of different algorithms at 2, 3, 4, and 5 levels.

Table 14 .
The fitness values obtained by different algorithms at higher threshold levels.

Table 15 .
Comparison of objective value, PSNR, SSIM, and FSIM values based on Tsallis entropy approaches.

Table 15 .
Comparison of objective value, PSNR, SSIM, and FSIM values based on Tsallis entropy approaches.