Adaptive Evolution of Metaheuristic Update Strategies Using Genetic Programming for Remote Sensing Image Fusion

1  Introduction

The swift progression of artificial intelligence (AI) technology has resulted in its growing implementation across various sectors, including healthcare, transportation, and industrial manufacturing. Notwithstanding this expansion, considerable hurdles persist, including technical impediments, data heterogeneity, and limitations dictated by particular application contexts. Comprehensive analyses of real-world problem models have revealed that these issues frequently display nonlinear, non-differentiable, and intricate properties. In this setting, meta-heuristic algorithms (MAs) have emerged as potent instruments for global optimization, providing a flexible and structure-agnostic approach to tackle these difficulties [1,2]. In contrast to conventional optimization methods, MAs are engineered to investigate intricate search spaces by emulating natural processes, biological evolution, or physical phenomena probabilistically. These algorithms utilize problem-specific knowledge to effectively pursue high-quality or near-optimal solutions. In recent years, meta-heuristic methodologies have exhibited exceptional efficacy across a diverse array of applications, encompassing wireless sensor networks (WSNs) [3,4], cloud computing resource allocation [5], communication security [6], intelligent manufacturing systems [7], machine learning [8], and numerous artificial intelligence endeavors [9,10].

The taxonomic classification of MAs is based on their foundational inspiration mechanisms, which are systematically divided into four principal categories as illustrated in Fig. 1. Metaheuristic algorithms are systematically categorized into four fundamental groups according to their sources of inspiration: Evolutionary-based algorithms, modeled on biological evolution principles, employ iterative population refinement through crossover, mutation, and selection mechanisms, with prominent implementations including Genetic Algorithm (GA) [11], Differential Evolution (DE) [12], Evolution Strategy (ES) [13], Biogeography-Based Optimization (BBO) [14], and Immunization Algorithm (IA) [15]; Swarm intelligence algorithms, derived from collective biological behaviors, leverage decentralized coordination strategies as exemplified by Particle Swarm Optimization (PSO) [16], Ant Colony Optimization (ACO) [17], Artificial Bee Colony (ABC) [18], Grey Wolf Optimization (GWO) [19], Whale Optimization Algorithm (WOA) [20], and Mountain Gazelle Optimization (MGO) [21]; Physics-based algorithms, grounded in physical/chemical laws, utilize force-directed solution space exploration techniques such as Simulated Annealing (SA) [22], Gravitational Search Algorithm (GSA) [23], Sine Cosine Algorithm (SCA) [24], and Multi-Verse Optimization (MVO) [25]; Human behavior-inspired algorithms, emulating cognitive-social dynamics, implement social learning paradigms including Teaching-Learning-Based Optimization (TLBO) [26], Harmony Search (HS) [27], Social Group Optimization (SGO) [28], and Special Forces Algorithm (SFA) [29]. This classification framework (Fig. 1) systematically organizes algorithmic design principles, wherein evolutionary methods formalize optimization through Mendelian genetics, swarm systems utilize stochastic agent-based models, physics-inspired techniques emulate natural forces, and human-centric approaches incorporate experiential adaptation mechanisms.

Figure 1: A brief classification of metaheuristic algorithms.

In recent years, research initiatives in MAs have focused on three primary trends. One of these paths entails the creation of wholly new algorithmic frameworks. The suggested Gannet Optimization Algorithm (GOA) [30], inspired by the foraging activities of gannets, employs U-shaped and V-shaped diving patterns for exploration and utilizes abrupt turns combined with random walks to enhance solutions. In a similar vein, Wang et al. [31] presented the Status-based Optimization (SBO) method, which emulates the human aspiration for status elevation. A different research direction focuses on the hybridization of various computer methods to collaboratively combine their respective benefits. Liang et al. [32] exemplified this approach by developing the hybrid IMPAEO algorithm, which integrates an enhanced marine predators algorithm with equilibrium optimization. Meanwhile, Song et al. [33] devised the HADEFP methodology through the strategic hybridization of differential evolution and flower pollination algorithms. A third research direction focuses on improving algorithmic efficacy by incorporating creative methodologies or parameter tweaking to boost performance across various problem domains [34].

Notwithstanding the formal diversity of metaheuristic algorithms, their fundamental evolutionary framework remains predominantly dependent on manually crafted update rules. The iterative formulations of algorithms, such as speed updates for Particle Swarm Optimization and position updates for Grey Wolf Optimization, are typically developed manually by researchers via empirical induction or phenomenological analogies. This technique demands a profound comprehension of the target problem’s characteristics and involves the adjustment of parameter weights following extensive trial-and-error experimentation. The constraints of these manual designs are becoming increasingly evident:

•   Empirical reliance and generalization deficiencies: Manual formulations are typically tailored for particular issue circumstances, and their sensitivity to parameters and inflexible strategies result in considerable performance variability when utilized across different domains. For instance, PSO’s fixed-rate update model often succumbs to local optimization in high-dimensional multi-peak problems, whereas GWO’s linear convergence mechanism finds it challenging to reconcile the demands of dynamic exploration and exploitation.

•   Absence of automation capability: The update protocols of current algorithms lack adaptive adjustment capabilities and are unable to dynamically generate or alter strategies during the optimization process. While hybrid algorithms can mitigate this problem, their foundational layers remain dependent on predetermined static rules, complicating the adaptation to abrupt modal shifts in intricate optimization landscapes.

•   Excessive utilization of computing resources: To enhance algorithmic robustness, researchers must regularly recalibrate formula parameters or devise novel hybrid solutions from the ground up. This requires considerable professional work and computational resources, hence limiting large-scale engineering applications.

To overcome these constraints, Meta-Black-Box Optimization (MetaBBO) has recently emerged as a potential approach, utilizing meta-learning to automate the design and adaptation of BBO algorithms with minimum expert involvement [35,36]. MetaBBO generally employs a bi-level framework, wherein a meta-level policy monitors the optimization dynamics of a low-level BBO process and adaptively formulates algorithmic designs based on performance feedback over a problem distribution. Current MetaBBO research may be classified as algorithm selection [37], algorithm configuration [37], solution manipulation [38], and algorithm generation [39]. Diverse learning paradigms, encompassing reinforcement learning [40], supervised learning [41], neuroevolution [42], and LLM-based in-context learning [43], have been investigated across various optimization contexts, including single-objective, multi-objective, multi-modal, large-scale, and multi-task optimization [44].

Notwithstanding the swift advancement of MetaBBO [45], the majority of current methodologies primarily depend on neural-network-based meta-policies—especially reinforcement learning [46]—to facilitate adaptive algorithm design via online state-action decision-making. While policy-learning-based frameworks exhibit robust performance, they generally necessitate meticulously crafted state representations, intricate reward functions, and substantial training interactions. This elevates algorithmic complexity and drastically diminishes interpretability. Furthermore, dependence on online policy inference incurs significant extra computing burden throughout the optimization process.

This article introduces a Genetic Programming-based Metaheuristic framework, referred to as GP-MAs, which emphasizes the evolving explicit symbolic update rules using genetic programming in an offline manner at the algorithmic level, in contrast to policy-learning-based MetaBBO frameworks. Genetic Programming (GP), a hyper-heuristic algorithm that autonomously generates computational programs [47], offers a potential foundation for tackling these challenges [48]. Rather of acquiring an implicit control policy, GP-MAs utilize the symbolic regression ability of GP to develop update formulas in an offline manner. The proposed system substitutes manually developed update rules with evolving symbolic expressions, therefore diminishing dependence on expert knowledge, preserving adaptability across various optimization settings, and reducing the online decision-making burden compared with neural MetaBBO approaches characteristic of neural MetaBBO models. In contrast to reinforcement-learning-based MetaBBO systems that primarily enhance implicit decision policies via online interaction, the proposed GP-MAs framework emphasizes the offline symbolic development of explicit update rules or search operators within metaheuristic algorithms. GP-MAs develop symbolic update rules utilizing predefined terminal sets (core MA variables) and function sets (mathematical/logical operators), assessing candidate formulations via a specified objective function to deliver adaptive update strategies for representative metaheuristic frameworks.

To implement and evaluate the efficacy of the proposed GP-MAs architecture, the Growth Optimizer (GO) [49], an algorithm derived from human learning and reflective behaviors, is utilized as a baseline human-inspired optimizer. This study selects GO as a sample human-inspired optimizer to assess the viability of the proposed GP-MAs system. This study enhances the GP-MAs framework by reconfiguring the learning-stage updating process of GO to create a hybrid variation called GPbasedGO. This integration incorporates GP-based symbolic strategy development into the learning phase of GO, establishing an adaptive symbolic updating mechanism. GPbasedGO produces a variety of candidate strategies using syntax tree encoding and employs a selection strategy considering convergence and exploration performance to equilibrate convergence velocity and exploratory capacity. To assess the efficacy of the proposed strategy, systematic tests are performed on the CEC2022 benchmark suite and a multispectral-panchromatic image fusion job. The benchmark tests assess convergence behavior and optimization performance, while the image fusion problem serves as an application case study of the proposed framework in typical remote sensing optimization scenarios.

The subsequent sections of this work are organized as follows. Section 2 thoroughly explains the fundamental approaches, including the Growth Optimizer (GO) algorithm, Genetic Programming (GP), and panchromatic-multispectral picture fusion techniques. Section 3 outlines the GP-MAs framework in conjunction with the architectural design and execution of the GPbasedGO algorithm. Section 4 delineates empirical findings and comparative studies. Ultimately, Section 5 consolidates the research findings and delineates potential avenues for future investigation.

2  Related work

This chapter will carefully outline the theoretical underpinnings and technical procedures relevant to this research, focusing on Growth Optimizer (GO), Genetic Programming (GP), and the key approaches to multispectral and panchromatic picture fusion. This material provides theoretical support for the following technique and experimental analysis.

2.1 Growth Optimizer

The GO algorithm, as a fledgling intelligent optimization approach, emulates the learning mechanism of reflection and self-improvement that is inherent in human development, with broad global search capabilities and considerable versatility. This section explains the core structure of the GO algorithm, including the initialization strategy, learning process, and the crucial notion of the reflection phase, which serves as a foundation for future algorithm modifications.

2.1.1 Initialization

The initialization process commences with the individuals of the 0th generation. For the ith individual Xi(0) in this generation, its jth dimension component xi,j(0) resides inside the interval [xi,jmin,xi,jmax]. Here, xi,jmin denotes the lower limit and xi,jmax signifies the upper bound for the jth dimension. The variable i denotes the ith individual in the population, while j signifies the jth dimension. The population size is represented by Np, while D signifies the number of dimensions. The population initialization is executed via a uniform distribution method, as illustrated in Eq. (1).

xi,j(0)=xi,jmin+rand(0,1)⋅(xi,jmax−xi,jmin)(1)

Eq. (1) demonstrates the initialization of each individual’s components within defined limits by a uniform random distribution, hence guaranteeing a varied beginning population across all dimensions and individuals.

2.1.2 Learning Phase

During the learning phase of the GO algorithm, individuals analyze their differences, explore the reasons behind these gaps, and learn from them, resulting in significant growth. To facilitate this process, four key types of gaps are defined: the difference between the leader and an elite individual (Gap1), between the leader and a lower-ranked individual (Gap2), between an elite and a lower-ranked individual (Gap3), and between two randomly selected individuals (Gap4). The mathematical model for these gaps is as follows:

{Gap1=x→best−x→betterGap2=x→best−x→worseGap3=x→better−x→worseGap4=x→L1−x→L2(2)

In this model, x→best denotes the leader of the group, while x→better represents an elite individual. Together, they form the upper tier of the social hierarchy. x→worse indicates a lower-ranked individual within the population. x→L1 and x→L2 are randomly chosen individuals, distinct from the i-th individual. Gapk (where k=1,2,3,4) signifies the difference between two individuals, offering learners valuable insights into these discrepancies for their benefit. It is important to note that within the current iteration (it) of the GO algorithm, to streamline the update and selection of individuals, GR should be sorted in ascending order.

The Learning Factor (LFk) is used to measure the strength of an individual’s learning from the kth set of gaps and is calculated as:

LFk=‖Gapk‖∑k=14‖Gapk‖,(k=1,2,3,4)(3)

where LFk reflects the relative size of the Euclidean distance of the group k gap among all gaps. The larger the LFk, the greater the impact of the gap on the individual’s learning, and the more the individual learns from that gap. For example, if Gap1 has a larger LF1, it indicates that the gap between leaders and elites is larger, and the individual needs to learn more from that gap to close the gap with leaders and elites.

The self-perception factor (SF) is used to assess an individual’s own situation and is calculated using the formula:

SF=GRGRmax(4)

where GR represents the individual’s growth resistance (with smaller values being better), and GRmax denotes the maximum growth resistance observed across all individuals. The scaling factor SF has values within the range [0,1]. The smaller the value of SF, the smaller the individual’s GR, indicating a higher position in the social hierarchy. Such individuals are more likely to engage in localized search (refined exploration mode). Conversely, the larger the value of SF, the greater the individual’s GR, placing them in a lower social tier. These individuals need to conduct more global searches (extensive exploration mode).

Knowledge Acquisition (KAk) is the core process in the learning phase, where individuals acquire knowledge from gaps through KAk, which is calculated as:

K→Ak=SF⋅LFk⋅Gap→k,(k=1,2,3,4)(5)

where KAk is the amount of knowledge an individual acquires from the kth set of gaps. SF and LFk together influence the size of KAk. For example, an individual in a higher social stratum (small SF) who learns more information (LFk) from a given gap will acquire a correspondingly larger amount of knowledge.

Individuals need to update their position after acquiring knowledge in the learning phase with the update formula:

x→iIt+1=x→iIt+K→A1+K→A2+K→A3+K→A4(6)

where xiIt is the position of individual i at the current iteration It, and xiIt+1 is the position of individual i at the next iteration It+1. Individuals grow by updating their position by absorbing the knowledge gained from the gap.

After updating, it needs to be verified whether the individual is actually progressing or not. If progress is made (i.e., the fitness value is lower), the updated individual is retained; if no progress is made, the new individual is retained with probability P2. This verification process ensures that individuals can truly improve themselves after learning and avoids the algorithm from falling into local optimality. The specific formula is shown in Eq. (7).

x→iit+1={x→iit+1if f(x→iit+1)<f(x→iit){x→iit+1if r1<P2x→iitelseelse(7)

where r1∈[0,1] is a uniformly distributed random number, and ind(i) denote the ranking of the ith individual based on an ascending order of the GR value. In cases where the individual fails to update, the retention of newly acquired knowledge is governed by the parameter P2, which is set to 0.001. Due to space limitations in Eq. (7), the full conditional statement for knowledge acquisition is expressed as: r1<P2 and ind(i)≠ind(1). This condition ensures that, even if an individual update fails, there is a small probability (0.001) that the learned information will be retained and the individual will be carried over to the next generation. Moreover, it guarantees that the current global best individual remains unchanged, thereby avoiding potential disruptions in the convergence of the algorithm.

Through the above process, the learning phase of the GO algorithm not only helps individuals to acquire knowledge from different gaps, but also adaptively adjusts the direction and intensity of learning through the self-awareness factor and the learning factor, so as to realize the balance between global search and local development.

2.1.3 Reflection Phase

The mathematical model of the reflection phase is described by Eqs. (8) and (9) as follows:

x→i,jit+1={lb+r4×(ub−lb)if r3<AFx→i,jit+r5×(Rj−x→i,jit)if r2<P3x→i,jitelse(8)

AF=0.01+0.99×(1−FEsMaxFEs)(9)

The reflection phase processes each dimension of the individual with three possible operations:

•   Keep the original value: If the random number r2≥P3, the individual’s jth dimension remains unchanged. This means that the individual has not learned anything new from other individuals in that dimension, or that the current state is considered good enough and does not need to be changed.

•   Updated under the guidance of a senior individual: If r2<P3, the individual’s j dimension is updated under the guidance of a senior individual, such as a leader or elite. The senior individual provides a goal direction toward which the individual moves, which helps the individual move closer to a more optimal solution in that dimension.

•   Reinitialization based on advanced individuals: If r3<AF, the jth dimension of the individual is reinitialized with a small probability. The process of reinitialization is to randomly select a new value in the search space, which increases the diversity of the population, prevents the algorithm from converging prematurely, and helps to jump out of local optima.

In Eq. (9), FEs is the current number of evaluations and MaxFEs is the maximum number of evaluations. As the algorithm iterates, AF gradually decays linearly from close to 1 to 0.01. At the beginning of the algorithm, AF is larger, and individuals have a higher probability of re-initialization, which helps in the global search. As the iteration progresses, AF gradually decreases, the probability of individual re-initialization decreases, and the algorithm gradually shifts to local search to refine the current solution.

The reflection phase helps people recognize and correct problems by motivating them to analyze and develop themselves. This stage encourages individual self-improvement so that the population as a whole can develop toward a better solution. Individuals can improve themselves by taking guidance from more experienced people.

GO was selected as a sample human-inspired optimizer due to its structured learning mechanism and strong performance in complex search contexts, making it a good baseline for evaluating the suggested methods. Although GO is utilized as the baseline optimizer in this study, the proposed framework is not limited to it and could be expanded to incorporate other population-based metaheuristics.

2.2 Genetic Programming

Genetic programming is a self-contained programming technique based on evolutionary ideas that creates program structures adapted to individual tasks by mimicking natural selection and genetic processes. In this study, genetic programming is used to produce key operators in optimization models. This section will present the individual generation method to genetic programming and its key evolutionary processes, laying the theoretical groundwork for the integration of GO algorithms with GP.

2.2.1 Individual Initialization

Individuals in genetic programming are usually encoded in a tree structure, where non-leaf nodes represent functions and leaf nodes are terminal symbols. A function can be found in the set of functions F={f1,f2,…,fn}, where fi includes operators (+,−,∗,/, Boolean operators, etc.) or functions (sin,tan,exp, etc.) as well as some expressions (cyclic expressions or conditional expressions, etc.). Terminators, on the other hand, can be selected from the set of terminators T={t1,t2,…,tm}, which can be variables or constants.

In the first stage of the algorithm, an initial population with diversity and a certain complexity must be generated, a process known as individual initialization. To generate an initial population with diverse structure, GP frequently uses the following initialization methods:

Full Method: This method uses function nodes to completely fill the non-leaf positions under the set maximum depth limit, and terminal nodes only appear in the leaf positions, generating an individual tree with regular structure and consistent depth.

Grow Method: This method randomly selects function or terminal nodes to fill in the node positions at each level, which makes the individual structure more irregular in depth and branching and is suitable for exploring the diversity of the structure space.

Ramped Half-and-Half Method: Combining the above two methods, alternating between the full-full method and the growth method to generate individuals in different depth ranges is the most commonly used initialization strategy, which is capable of balancing the structural diversity of the population with the search breadth.

2.2.2 Basic Operations of Genetic Programming

The main evolutionary process of genetic programming is made up of three basic operations: selection, crossover, and mutation, which mimic the evolutionary process of organisms and drive the ongoing optimization of programmed people under the guidance of the fitness assessment mechanism.

(1) The selection operation selects which people have the opportunity to enter the next generation and reflects the supremacy of genetic programming. Its primary purpose is to identify highly adaptable programmed individuals based on their particular fitness levels. Commonly used selection strategies are:

Roulette Wheel Selection: The probability of an individual being selected is directly proportional to its fitness value, and is suitable for use in situations where there is a large difference in fitness.

Tournament Selection: A number of individuals are randomly selected from the population for a tournament to select the one with the highest fitness, which has a strong selection pressure.

Rank Selection: Assigning probabilities according to the order of fitness to avoid monopolization of the breeding process by over-adapted individuals.

Adjusting the selection strategy appropriately can avoid the rapid loss of population diversity while maintaining the inheritance of superior genes.

(2) Crossover operation: Crossover is the main ?recombination? mechanism in genetic programming, which aims to combine the structural information of two good individuals to create potential high-quality offspring. As illustrated in Fig. 2. The basic process is as follows:

•   Randomly select a subtree (substructure) from each of the two parents.

•   Swap these two subtrees.

•   Generating two new individuals as offspring.

Figure 2: Sample diagram of crossover operation.

(3) Mutation operation: The mutation operation introduces a small structural disturbance to prevent the algorithm from falling into the local optimum and maintain the diversity of the population. As illustrated in Fig. 3, in genetic programming, the most commonly used is subtree mutation:

•   Randomly select a subtree in the individual as the mutation point.

•   Replace the position with a newly generated random subtree.

•   Form a new program structure.

Figure 3: Sample diagram of mutation operation

In addition, there is also point mutation, which is to randomly replace a node in the individual (such as replacing a function symbol or variable symbol), which is less destructive but introduces less information. The mutation probability is usually set small to avoid destroying the existing excellent structure.

2.3 Panchromatic and Multi-Spectral Image Fusion

In the realm of remote sensing image processing, multi-source data fusion is an important technique for increasing image quality. Among these techniques, the fusion of panchromatic (Pan) and multispectral (MS) pictures is especially important. Fig. 4 shows a schematic diagram of the pan-sharpening process. By combining multispectral images with high spectral but low spatial resolution and panchromatic images with low spectral but high spatial resolution, the fusion process creates an image with both high spectral and spatial resolution. Remote sensing image fusion technology encompasses a variety of established methods, typically categorized as follows:

1.   Color Transformation Methods: The Intensity-Hue-Saturation (IHS) transformation represents a widely adopted technique. By converting RGB color space into IHS space, this method replaces the intensity component with the high-resolution panchromatic image, thereby enhancing spatial resolution while preserving spectral content.

2.   Statistical Analysis Methods: Principal Component Analysis (PCA) is a representative statistical method that decorrelates spectral bands into orthogonal principal components. The first principal component, which captures the largest variance, is typically replaced by the panchromatic image. While effective in maintaining spectral characteristics, PCA may fall short in capturing spatial details.

3.   Numerical Methods: The Brovey transformation exemplifies a numerical approach, leveraging arithmetic operations within the RGB space to inject spatial information from the panchromatic image. It modulates the brightness component of the multispectral image while retaining its spectral fidelity.

4.   Hybrid Methods: Hybrid fusion techniques integrate multiple fusion paradigms to leverage their respective strengths. The proposed fusion strategy in this study is a representative hybrid approach, integrating intelligent optimization algorithms with conventional fusion techniques to enhance overall performance. In this study, image fusion is used as an optimization-driven application scenario to evaluate the effectiveness of the proposed metaheuristic learning framework rather than proposing a new fusion algorithm.

Figure 4: Schematic diagram of panchromatic image and multispectral image fusion.

Image fusion methods based on parameter optimization have emerged as a significant research focus. The core idea lies in constructing a mathematical mapping between the panchromatic and multispectral images. The general model is described as:

F=M+g(P−D)(10)

where F is the fused high-resolution image, M denotes the multispectral band, P is the panchromatic image, D represents the detail component, and g is the mixing parameter that controls the balance of spatial and spectral contributions.

In the IHS framework, the intensity component I is formulated as:

I=∑i=1nαiMi(11)

Here, αi is a non-negative coefficient. In RGB images, αi is commonly set to 1/3 to ensure equal contribution from each channel. For multispectral images with n bands, αi is generally set to 1/n to maintain balanced spectral contribution across all channels.

To further reduce spectral distortion, the goal is to ensure that the panchromatic image closely approximates a linear combination of the multispectral bands:

P≈∑i=1nαiMi(12)

This leads to an optimization problem that minimizes the spectral discrepancy:

|P−∑i=1nαiMi|2(13)

The optimal αi can be determined via the Least Squares Method (LSM):

α=(MTM)−1MTP(14)

This optimization framework effectively enhances the spatial resolution of multispectral images while preserving their spectral integrity, thereby supporting high-quality remote sensing applications.

In image processing, quantitative evaluation of fusion performance is crucial for verifying algorithm efficacy. Multiple evaluation metrics are commonly employed to assess spatial enhancement, spectral preservation, and visual fidelity.

1. Root Mean Square Error (RMSE): RMSE is a fundamental metric that quantifies pixel-level deviations between the fused image and the reference. It is defined as:

RMSE=1L×W∑i=1L∑j=1W(Fi,j−Mi,j)2(15)

where L and W denote the image dimensions, F and M represent pixel intensities of the fused and reference images, respectively. Lower RMSE indicates better fidelity to the reference image.

2. Relative Dimensionless Global Error in Synthesis (ERGAS): ERGAS evaluates the relative discrepancy between fused and reference multispectral images, accounting for resolution differences:

ERGAS=100×1n∑i=1n(RMSEiμi)2(16)

where n is the number of bands and μi is the mean of the ith band. Lower ERGAS implies higher overall image quality.

3. Spectral Angle Mapper (SAM): SAM measures spectral similarity using vector angles:

SAM=arccos⁡(∑iMi⋅Fi∑iMi2⋅∑iFi2)(17)

Smaller SAM values reflect better spectral preservation.

4. Relative Average Spectral Error (RASE): RASE quantifies the mean spectral deviation relative to each band:

RASE=100μM1n∑i=1nRMSEi2(18)

Lower RASE values indicate superior spectral consistency.

5. Universal Image Quality Index (UIQI): UIQI integrates brightness, contrast, and structural similarity:

UIQI=4σMF(σM2+σF2)⋅2μMμF(μM2+μF2)(19)

UIQI values range from −1 to 1, with values closer to 1 indicating better overall quality.

6. Correlation Coefficient (CC): CC assesses linear correlation between the fused and reference images:

CC=∑i,j(Mi,j−μM)(Fi,j−μF)∑i,j(Mi,j−μM)2∑i,j(Fi,j−μF)2(20)

CC values close to 1 suggest strong similarity and effective information preservation.

Collectively, these evaluation metrics provide a robust foundation for comprehensive performance assessment in remote sensing image fusion, facilitating the development of more accurate and reliable fusion methodologies.

3  GP-MAs Framework and GO Algorithm Based on the GP-MAs Framework

This chapter provides a detailed introduction to the GP-MAs framework. It begins by outlining the framework’s overall architecture and then goes into detail into the design of the relevant parameters. It also investigates the integration of Genetic Programming-based Metaheuristic Algorithms (GP-MAs) into the Growth Optimizer (GO) and gives a thorough examination of the updating mechanisms specifically designed for GO.

3.1 Overall Structure

Fig. 5 shows the architecture of the proposed GP-MAs-GO system, with red lines indicating data flow. This study’s basic algorithm is the original Growth Optimizer (GO). The GP-MAs module is only included in the learning phase of GO, with the goal of constantly improving the learning method. The reflection phase and overall structural integrity of GO remain unaffected. In Step 1, researchers manually define the function set, terminal set, and objective function for evaluating the established learning algorithms. These capabilities are subsequently incorporated into the GP-MA system. Within the GP-MAs framework, a population of genetic programming (GP) individuals is established, with each representing a potential learning formula. Advantageous learning mechanisms are preserved by evolutionary processes like as selection, crossover, and mutation, while less successful ones are gradually destroyed.

Figure 5: System framework diagram.

The assessment process for each individual in the GP population comprises putting the candidate formula into the learning phase of GO and replacing the original update rule. The improved GO method is then applied to a set of benchmark functions, with the ultimate optimization performance serving as the fitness value for each individual. This ensures that only learning strategies that improve GO are kept.

In Step 2, the original handcrafted formula in GO is replaced with the most effective evolved formula obtained using the GP-MAs approach. The resulting algorithm, GPbasedGO, distinguishes itself from classic GO by using a dynamically changeable learning mechanism. In contrast to the fixed, human-engineered updating rules used in traditional GO, GPbasedGO adapts its learning method to changing issue landscapes and evolutionary stages.

The growth of learning approaches via GP-MAs takes place offline, allowing the best formula to be used across several optimization issues without retraining. In practical optimization research, the manual creation of evolutionary formulas typically requires iterative parameter adjustments, variable alignments, heuristic experimentation, and comprehensive performance evaluations, resulting in significant instability and difficult cost assessment. The proposed GP-MAs framework replaces human design patterns with automated symbolic evolution, with the GP process performing variable matching and performance assessment on a continuous basis. While the offline GP evolution phase may require additional computational time, the created symbolic methodologies can be reused in future optimization projects, reducing redundant manual redesign efforts and extended human participation. Thus, the proposed framework provides a more adaptable and automated approach for symbolic strategy building in metaheuristic optimization.

Unlike traditional GP-based hyper-heuristics, which evolve entire algorithms or heuristics, GP-MAs evolve symbolic update rules nested within a fixed metaheuristic backbone. This constraint minimizes search space complexity while preserving the interpretability of the generated operators. In contrast to MetaBBO approaches, which train policy networks or selection strategies, GP-MAs directly optimize explicit mathematical update expressions.

3.2 The GP-MAs Framework

The GP-MAs (Genetic Programming-based Metaheuristic Algorithms) framework employs genetic programming to automatically generate update techniques for metaheuristic optimization algorithms. The system encodes search behaviors in adaptive mathematical expressions, which are then iteratively refined to increase optimization performance. Table 1 provides an overview of GP-MA parameter configurations. The GP parameter settings used in this study were based on commonly acknowledged configurations in the GP literature and preliminary empirical findings. These parameter sets intended to maintain an equilibrium between evolutionary diversity and computational stability as symbolic strategies evolved. Nonetheless, the efficiency of GP-based symbolic evolution may be influenced by hyperparameter selection. A more extensive sensitivity analysis of Gaussian Process hyperparameters will be addressed in future studies.

The detailed steps of the GP-MAs framework are as follows:

1.   Initialization: A population of individuals is randomly initialized using a predefined set of functions and terminals. To enhance structural diversity, half of the individuals are generated using a depth-first growth method, while the remaining half use breadth-first construction.

2.   Expression Construction and Evaluation: Each individual encodes a candidate metaheuristic component (e.g., an update rule or control policy) as a tree structure. These individuals are converted into executable expressions via in-order traversal and then embedded into a target metaheuristic algorithm, where their performance is evaluated using a predefined objective function.

3.   Elitism: The top-performing individuals, based on fitness scores, are preserved in the next generation to ensure the retention of high-quality solutions.

4.   Variation Operators: The remaining individuals undergo variation through crossover and mutation, with parent selection guided by a roulette-wheel selection mechanism. Structural constraints are applied to maintain the syntactic and semantic validity of the expressions.

5.   Fitness Evaluation: The fitness of the new individuals is recalculated by integrating them into the algorithm and re-evaluating their performance.

6.   Termination: Steps 2 to 5 are repeated until the stopping criterion is met, such as reaching the maximum number of generations or achieving a satisfactory performance level.

Algorithm 1 outlines the workflow of the GP-MAs framework. Inspired by prior studies on evolutionary formula generation in metaheuristics, this approach ensures structural diversity at initialization via a balanced tree construction strategy. Each individual, represented as a candidate metaheuristic rule, is transformed into an executable formula and evaluated within a host algorithm. Top individuals are retained via elitism, while the remaining ones undergo crossover and mutation under structural constraints. The evolution proceeds iteratively until convergence or reaching the iteration limit.

To avoid excessive growth of GP trees (bloat), the maximum tree depth was restricted during the evolutionary process according to the predefined depth constraints in Table 1. Structural constraints were also applied during crossover and mutation to maintain syntactic validity and control expression complexity. Since the GP evolution stage is performed offline, the additional computational overhead is incurred only once during strategy generation and does not affect the subsequent online optimization stage. The overall computational complexity of GP-MAs can be approximated as O(GGP×PGP×Nf×R×FEs×Cf), where GGP denotes the number of GP generations, PGP is the population size, Nf represents the number of training functions, R is the number of independent runs, and Cf denotes the computational cost of fitness evaluation. In addition, protected division and bounded nonlinear operations were adopted during GP evolution to improve numerical stability and reduce the risk of singularities or excessively large update values caused by nested trigonometric operators.

To quantitatively assess the performance of evolved formulas, a customized evaluation function F is designed for the GP-MAs framework. To provide representative optimization landscapes during GP evolution, four benchmark functions are selected from the CEC2022 benchmark suite: F1 (unimodal), F3 (basic multimodal), F7 (hybrid), and F10 (composition). These functions were selected to represent different optimization characteristics, including unimodal, multimodal, hybrid, and composition landscapes, thereby improving the diversity of training scenarios during GP evolution. Each evolved formula is integrated into a metaheuristic algorithm variant, referred to as new-MAs, and tested on these functions.

Let fitj be the mean result of new-MAs on function j over five independent runs, and fitj∗ denote the known global optimum. The overall fitness score is computed as:

Fitness=∑j∈{1,3,7,10}fitjfitj∗(21)

This aggregated metric provides an indirect assessment of the transferability and adaptability of evolved strategies across representative optimization scenarios.

Although the chosen benchmark functions include representative unimodal, multimodal, hybrid, and composition landscapes, the current study does not address extremely dynamic, limited, or non-stationary real-world optimization contexts. As a result, the current studies should be interpreted as representative validation rather than universal generalization verification. Future research will focus on more difficult optimization scenarios. The computational complexity of GP-MAs is dominated by GP evolution, which runs at O(G × P × F), where G is the number of generations, P is population size, and F is the fitness evaluation cost. Because GP evolution is performed offline, its cost has no impact on the runtime complexity of GPbasedGO during optimization.

3.3 The GP-GO Algorithm

To minimize numerical instability and overly complicated structures, evolving expressions are controlled during GP search with protected operators, bounded depth control, and fitness-based pruning. Only expressions that meet predefined complexity constraints are kept. This study proposes GPbasedGO, an upgraded form of the Growth Optimizer (GO) algorithm optimized using the GP-MAs framework. Within this paradigm, GP-MAs independently evolve the learning rule. The evaluation function F is set to a size of 10. In accordance with recognized techniques, four sample functions (F1, F3, F7, and F10) from the CEC2022 benchmark are chosen for evaluation.

The terminal set 𝒯 and function set ℱ critically determine GP performance. For GO optimization, these sets are designed as:

{𝒯={Gap1,Gap2,Gap3,Gap4,SF}ℱ={+,−,×,÷,norm,sin,cos,tan}

To evaluate the contribution of GP-evolved operators, several ablation variants are defined: GO-original: baseline Growth Optimizer; GO-SF: only self-perception factor enabled; GO-fixed-rule: handcrafted update rule; GPbasedGO: full proposed model;

𝒯 encodes social hierarchy dynamics: Gap1−4 quantify rank disparities (e.g., leader-elite, elite-lower, random pairs) to guide knowledge acquisition through hierarchical learning mechanisms. SF (self-perception factor) dynamically scales learning intensity via growth resistance (GR), balancing local refinement (SF≪1) and global exploration (SF≫1). This adaptive scaling mirrors parameter self-adaptation strategies in differential evolution.

ℱ integrates arithmetic operators for evolutionary recombination and norm for numerical stability. Trigonometric functions inject nonlinear perturbations to escape local optima while maintaining diversity. This hybrid design balances deterministic hierarchical learning with stochastic exploration, accelerating cnvergence while preserving solution diversity. Embedded self-adaptive components reduce parameter sensitivity while maintaining evolutionary pressure toward optimal regions, aligning with the “learning-reflection-self-improvement” philosophy of human-inspired optimization.

The update rule for this process can be expressed as follows:

Term1=‖Gap1‖⋅Gap1+‖Gap2‖⋅Gap2,Term2=‖Gap3‖⋅Gap3+‖(Gap1−Gap1)−(SF⋅SF)‖⋅Gap4

The numerator of the update formula can then be written as:

Numerator=SF⋅(Term1+Term2).(22)

The denominator is given by:

Denominator1=‖cos⁡(tan⁡(cos⁡(‖Gap1‖)))‖+‖Gap2‖,Denominator2=‖Gap3‖+‖sin⁡(‖cos⁡(Gap3−Gap2)‖)‖,(23)

Thus, the complete update equation becomes:

Update=NumeratorDenominator1+Denominator2.(24)

Compared with the original formulation, the new update equation is automatically evolved through Genetic Programming (GP), enabling more adaptive information interaction among different gap terms. Specifically, Gap1 and Gap2 mainly emphasize elite-guided exploitation, while Gap4 introduces stochastic diversity to enhance global exploration. The self-perception factor SF adaptively adjusts the learning intensity according to the evolutionary state of individuals. In addition, the nonlinear operators, including trigonometric functions and normalization terms, increase directional diversity and help maintain the balance between exploration and exploitation during the search process. The denominator terms further normalize the update magnitude, which contributes to numerical stability and avoids excessively large search steps. Therefore, compared with the original manually designed formulation, the evolved symbolic update rule provides a partially interpretable adaptive learning mechanism for complex optimization tasks.

3.4 Fusion Method Based on GPbasedGO Algorithm

To further improve the performance of image fusion, the basic fusion formula (Eq. (10)) is modified, and the updated formulation is expressed in Eq. (25):

Fi=Mi+giΓi(P−I)(25)

where Γi denotes a weight matrix that plays a pivotal role in determining the relative contributions of the panchromatic image P and the low-resolution multispectral image Mi during the fusion process, thereby directly influencing the quality of the fused output. The definition of Γi is given in Eq. (26):

Γi=β(i)ΓP+(1−β(i))ΓM(26)

where β(i) is a coefficient associated with the i-th spectral band of the multispectral image. By adjusting β(i), the relative proportions of ΓP and ΓM within Γi can be flexibly controlled. The computation of ΓP and ΓM is described in Eqs. (27) and (28), respectively:

ΓP=exp⁡(−λ‖∇P‖4+ε)(27)

ΓM=exp⁡(−λ‖∇Mi‖4+ε)(28)

where ∇P and ∇Mi denote the gradients of the panchromatic and multispectral images, respectively. The gradient encodes pixel intensity variations, enabling the preservation of edge and detail features during fusion. Incorporating gradient information thus facilitates more precise processing in the fusion stage. The term ε is a small constant introduced to prevent division by zero, ensuring numerical stability and accuracy. The denominators ‖∇P‖4+ε and ‖∇Mi‖4+ε are safeguarded against zero values. The parameter λ governs the degree of smoothness in the resulting fused image. Empirical evaluation demonstrates that setting λ=10−9 and ε=10−10 yields substantial improvements in spectral resolution. Under this configuration, the fused image more faithfully preserves spectral characteristics across different bands, thus enhancing spectral fidelity and meeting the high-quality imaging requirements of remote sensing and related applications.

Through the proposed definition of Γi and careful parameter tuning, the integration of multispectral and panchromatic imagery can be performed more effectively, producing high-quality fusion results. This provides superior input data for subsequent remote sensing analysis and processing tasks. The overall workflow of the proposed GP-GO-based remote sensing image fusion approach is illustrated in Fig. 6.

Figure 6: Flowchart of the proposed fusion method

4  Experimental Results

This section describes the experimental framework for testing the proposed GP-based GO algorithm. Representative benchmark functions from the CEC2022 test suite are used to evaluate optimization performance in difficult settings. A remote sensing image fusion task (panchromatic-multispectral fusion) further validates the practical efficacy of GP-basedGO. This dual-assessment framework analyzes robustness, adaptability, and optimization skills in both synthetic and real-world environments.

4.1 Experimental Setup

A comprehensive comparative study is conducted against fifteen algorithms, comprising classical methods (PSO, GA, ASO, BOA, SA, CS, GWO, DE, CRO) and recent state-of-the-art metaheuristics (GO, WOO [50], EMWPA [51], FNO [52], CHO [53], LEE [54]). To maintain fairness, all algorithms use the same population sizes and fixed fitness evaluations (FES). Each algorithm is run 30 times individually on each benchmark function to mitigate the effect of initialization randomization on solution accuracy and assure statistical reliability. The performance is evaluated using the mean and standard deviation of the answers, as well as the Wilcoxon rank-sum test. All experiments in this study were carried out on a Windows 11 computer platform outfitted with an NVIDIA GeForce GTX 1650 GPU with 4 GB of video RAM and CUDA parallel computing capabilities. The central processing unit (CPU) is an Intel Core i7-9750H processor with six cores running at 2.60 GHz, and the system comes with 16 GB of RAM.

To assess the performance of the proposed GPbasedGO algorithm, the standard CEC2022 benchmark suite was used. This test set includes four types of functions: unimodal, basic multimodal, hybrid, and composition functions. The unimodal function (F1) is specifically developed to evaluate the algorithm’s local search capabilities and convergence accuracy because it only has one global optimum.

The basic multimodal functions (F2–F5) contain a significant number of local optima, which are used to assess the algorithm’s global exploration capacity and ability to avoid premature convergence to local minima. These roles present a considerable challenge in terms of striking an efficient balance between exploration and exploitation.

The hybrid functions (F6–F8) are created by merging various basic functions with different properties, resulting in more complex and deceiving landscapes. These functions test the algorithm’s robustness and flexibility in dynamic and heterogeneous optimization contexts.

Finally, the composition functions (F9–F12) are the most complex category. They are created by combining many hybrid or basic functions with different biases, rotations, and scaling transformations. These functions add considerable levels of nonlinearity, modality, and search space irregularity, serving as a demanding testbed for evaluating the algorithm’s capacity to locate and transition between numerous diverse regions of the search space.

Through extensive experimentation across this broad range of benchmark functions, the GPbasedGO algorithm’s effectiveness and generalization potential may be comprehensively tested in optimization situations of different difficulty and landscape complexity.

4.1.1 Analysis of 10-Dimensional Results

Tables 2–4 compare the GP-basedGO method to fifteen different optimization strategies for 10-dimensional problems. The comparison of mean values demonstrates that GP-basedGO has greater overall performance. The statistics demonstrate that GP-basedGO generally outperforms other algorithms, with its mean metric exhibiting statistically significant improvements on a wide range of test functions, particularly F1, F4, F6, and F11. However, for some functions, like as F3 and F7, a few rival algorithms may perform similarly to or slightly better than GPbasedGO.

When GPbasedGO is compared to classical algorithms (such as PSO, GA, and DE), its distinct advantages become clear. According to the tables, GP-basedGO outperforms other methods in a variety of test conditions. When compared to PSO, GP-basedGO has smaller mean values and faster convergence across functions, indicating improved performance in optimizing high-dimensional complex landscapes. When compared to GA, GP-basedGO has a substantially lower standard deviation, indicating improved stability, precision, and reduced sensitivity to initial circumstances and random causes. GP-basedGO outperforms DE and related algorithms in terms of multimodal function search. Furthermore, its advantages over algorithms such as GWO and CS are essentially based on balancing unimodal and multimodal optimization, demonstrating its strong problem-solving abilities for complicated problems.

Furthermore, as compared to current state-of-the-art metaheuristics, GP-basedGO maintains a significant competitive advantage. In comparison to high-performing approaches such as CHO and EMWPA, GP-basedGO outperforms them on complicated multimodal landscapes, where it more successfully avoids local optima to ensure higher precision. While recently presented algorithms such as FNO and LEE excel at certain, isolated functions, GP-basedGO provides a far more balanced and versatile performance profile over the whole benchmark suite, demonstrating its broad application.

Fig. 7 shows the convergence curves of various methods in a 10-dimensional space. The picture depicts representative portions of convergence behavior, with subplot names identifying the relevant benchmark functions. It is clear that GP-basedGO (shown by the red curve) has strong optimization capabilities, with a constant downward trajectory that eventually outperforms practically all other algorithms. Most algorithms converge quickly within the first 500 evaluations; however, on complex functions like F2 and F12, competitors such as CS and CRO tend to plateau early on. However, GP-basedGO well balances local exploitation and global exploration, which is an important property for optimizing complex systems. This equilibrium avoids premature convergence to local optima while efficiently traversing the whole search space to find high-quality solutions.

Figure 7: 10-dimensional convergence plots of the algorithm.

4.1.2 Analysis of 20-Dimensional Results

Tables 5–7 compare the 20-dimensional outcomes of the GP-basedGO method to fifteen other optimization algorithms. Notably, the GP-based GO updating approach was developed with only 10-dimensional training data. Its higher performance on 20-dimensional test functions reveals high dimensionality transferability, implying that the evolved formula reflects scalable optimization processes rather than simply overfitting the low-dimensional training set. The average values in the tables show that GPbasedGO outperforms other methods on the majority of functions. On specific test functions such as F1, F2, F8, and F12, GP-basedGO outperforms other algorithms in terms of average value metrics. However, for some functions, such as F10, a few rival algorithms may perform similarly or slightly better than GPbasedGO.

When compared to standard algorithms (e.g., PSO, GA, DE, GWO, and GO), GP-basedGO has distinct advantages and beats them in stability. When compared to PSO, GPbasedGO has lower average values on most functions and faster convergence in the initial stage, indicating improved performance in high-dimensional complex function optimization. In comparison to GA, GP-basedGO has a lower standard deviation, indicating more stable results with less effect from initial conditions and random causes. GP-basedGO outperforms DE in terms of multimodal function search. Furthermore, it beats algorithms like as GWO and GO in both unimodal and multimodal settings. Furthermore, GP-basedGO remains very competitive with newer cutting-edge metaheuristics (e.g., WOO, FNO), indicating its robust capacity to tackle complicated problems in a variety of terrains.

Fig. 8 shows the convergence curves of various methods in a 20-dimensional space. We offer only representative parts of the convergence graphs here, with subplot names specifying the corresponding benchmark functions. As illustrated by the convergence curves, the move to 20 dimensions causes many algorithms to plateau due to the “curse of dimensionality” (e.g., CS and CRO on specific functions). However, GPbasedGO has a consistent decreasing tendency, displaying significant optimization abilities and surpassing almost all other algorithms in the end results. This robustness comes from the GP-evolved formula’s inherent scalability, which uses non-linear operators to generate various step sizes, eliminating diversity loss and assuring robust performance in higher-dimensional spaces. GPbasedGO iterates between two phases, each having a growth cycle, allowing it to efficiently balance local and global exploration. During the seedling growth stage, the population update approach slows and stabilizes the search space exploration process. This strategy allows for some exploration of new regions while focusing more on the in-depth evolution of the surrounding environment, preventing premature convergence to local optima and enabling extremely effective search space utilization.

Figure 8: 20-dimensional convergence plots of the algorithm.

4.2 Experimental Evaluation of GPbasedGO on Multispectral and Panchromatic Image Fusion

This section looks at the performance and applicability of the Genetic Programming-based Growth Optimizer method (GPbasedGO) for adaptive multispectral and panchromatic image fusion applications. To test the usefulness of GP-basedGO in tackling such fusion difficulties, benchmark experimental data from two satellites with different theme characteristics are used. The publicly available NBUPansharpRSData dataset is used, which can be found at https://github.com/starboot/NBUPansharpRSData. This comprehensive dataset contains multispectral photos acquired by satellites including QuickBird, IKONOS, GF, and WorldView.

Four sample scenarios from the popular QuickBird and GF satellite pictures are chosen for focused and representative examination. The experimental procedure is organized in the following five topic modules:

1.   Dataset and Parameter Configuration: Description of the remote sensing data characteristics and algorithm parameter settings.

2.   Urban Scene Analysis: Performance comparison of 15 state-of-the-art fusion algorithms in urban landscape scenarios.

3.   Vegetation Scene Analysis: Quantitative evaluation of 15 benchmark algorithms in green vegetation areas.

4.   Aquatic Environment Analysis: Objective assessment of 15 established methods in water-dominated environments.

5.   Coastal Region Analysis: Comprehensive comparison of 15 representative algorithms in coastal terrain conditions.

4.2.1 Experiments in Urban Scenarios

Table 8 shows that the suggested GP-based GO consistently beats all four urban settings. In QuickBird (Image 1), GPbasedGO gets the best values in RASE (20.39585), RMSE (12.63616), UIQI (0.976627), and CC (0.970901), while maintaining an ERGAS value of 382.4032, which is extremely close to the optimal result attained by EMWPA. Although EMWPA outperforms GP-basedGO in ERGAS and SAM, the proposed technique has greater spatial-spectral consistency, as evidenced by higher UIQI and CC values.

In the challenging QuickBird (Image 2) scene, GP-basedGO beats all other approaches by achieving the best results on all six evaluation metrics: ERGAS, SAM, RASE, RMSE, UIQI, and CC. In particular, when compared to the original GO, the proposed method reduces ERGAS from 27.77494 to 25.39373 while increasing CC from 0.752443 to 0.828362, indicating that the introduced Gaussian-process-guided search strategy significantly improves the ability to preserve structural details in complex urban textures.

In both the GF (Image 1) and GF (Image 2) datasets, GPbasedGO consistently outperforms other approaches. GF (Image 2) outperforms all baselines, with significant decreases in ERGAS, RASE, and RMSE, as well as the highest UIQI and CC values. These results show that the proposed technique may effectively balance spectral fidelity and spatial detail preservation despite differences in sensor qualities and scene complexity.

Visual comparisons in Fig. 9 demonstrate the efficacy of GP-based GO. The proposed method yields fused images with improved edge structures, richer building textures, and more natural color distributions, especially at road crossings, rooftops, and heavily populated regions. Compared to the low-resolution multispectral and panchromatic inputs, the fusion outputs have higher contrast and better spatial characteristics without causing significant spectral distortion. This observation supports the better quantitative scores in Table 8, indicating that GP-based GO provides robust and reliable fusion performance for urban remote sensing scenarios.

Figure 9: Visual assessment for marked region in Urban dataset.

4.2.2 Experiments in Green Vegetation Scenarios

Table 9 shows that the proposed GP-basedGO outperforms competing approaches in green vegetation conditions across four test cases. For QuickBird (Image 1), GPbasedGO achieves the best values in five out of six evaluation measures, including ERGAS (15.68661), RASE (53.17873), RMSE (36.47474), UIQI (0.849425), and CC (0.985096), with a SAM value that is just somewhat inferior than that of EMWPA. This shows that the suggested technique is extremely effective in preserving vegetative structures while maintaining global spectral uniformity.

In the tough QuickBird (Image 2) scene, GPbasedGO performs optimally on all six quality measures. Compared to the original GO, the proposed technique lowers ERGAS from 16.60697 to 16.56209 and SAM from 14.60079 to 14.44333, while increasing UIQI and CC to 0.866672 and 0.954173, respectively. These enhancements demonstrate that the Gaussian-process-guided search mechanism improves the optimizer’s capacity to capture fine-grained vegetation textures, edge transitions, and small spectral fluctuations in complex green regions.

The GF (Image 1) and GF (Image 2) datasets show a similar trend, with GPbasedGO achieving the best overall performance. In particular, for GF (Image 2), the proposed technique consistently achieves the lowest distortion-related metrics and the highest similarity-related metrics among all analyzed algorithms, indicating strong robustness under a variety of sensor and vegetation distributions. This higher performance indicates that GP-basedGO may efficiently balance spatial detail augmentation and spectral fidelity preservation in heavily textured vegetation images.

Fig. 10 shows visual comparisons that confirm the quantitative conclusions. GP-basedGO fusion findings show crisper canopy borders, sharper leaf-cluster textures, and more natural color transitions in densely vegetated areas. Furthermore, the suggested method preserves local contrast in shaded regions and complicated background boundaries, thereby avoiding the texture blurring and spectral inconsistencies seen in traditional optimization-based fusion methods. Visual observations are highly consistent with superior objective scores in Table 9, thus demonstrating the usefulness and robustness of GP-basedGO for remote sensing image fusion in green vegetation contexts.

Figure 10: Visual assessment for marked region in green vegetation dataset.

4.2.3 Experiments in Water Scenarios

The quantitative results summarized in Table 10 further demonstrate the superior fusion capability of the proposed GPbasedGO in water scenarios. Across the four test cases, GPbasedGO consistently achieves the best or highly competitive performance, confirming its robustness in handling homogeneous regions and subtle spectral variations that are particularly challenging in water-dominated scenes.

For QuickBird (Image 1), GPbasedGO achieves the best values in SAM, RASE, RMSE, UIQI, and CC, while its ERGAS value remains highly competitive and only marginally inferior to GO. In particular, the proposed method significantly improves SAM to 9.0233 and CC to 0.8870, indicating that GPbasedGO is highly effective in preserving spectral smoothness and boundary consistency in water-land transition regions.

In the QuickBird (Image 2) scene, GPbasedGO delivers the optimal performance across all six evaluation metrics, clearly outperforming all competing methods. Compared with the original GO, the proposed method reduces ERGAS from 33.0738 to 33.0507 and improves CC from 0.9583 to 0.9587. These improvements, although numerically moderate, are highly meaningful in water scenarios where local textures are relatively weak and spectral distortions are more likely to dominate the fusion quality. This confirms that the Gaussian-process-guided search mechanism can more effectively capture subtle intensity gradients and preserve fine shoreline structures.

A similar superiority can be observed in the GF (Image 1) and GF (Image 2) datasets. In GF (Image 1), GPbasedGO achieves the best performance on all six metrics, demonstrating excellent detail preservation in smooth water surfaces and surrounding land boundaries. For GF (Image 2), the proposed method again attains the best values in ERGAS, RASE, RMSE, UIQI, and CC, while its SAM remains highly competitive and only slightly inferior to LEE. These results highlight the strong generalization capability of GPbasedGO under different sensor characteristics and diverse water distributions.

The visual comparisons in Fig. 11 further support the quantitative analysis. The fusion results generated by GPbasedGO exhibit clearer shoreline contours, more natural tonal transitions over water surfaces, and better preservation of weak boundary structures in shadowed or reflective regions. Compared with the baseline methods, the proposed approach effectively suppresses spectral distortion and avoids the over-smoothing artifacts that commonly occur in homogeneous water areas. The high visual consistency between water bodies and adjacent land regions is in strong agreement with the superior UIQI and CC values reported in Table 10, further validating the effectiveness of GPbasedGO for remote sensing image fusion in water scenarios.

Figure 11: Visual assessment for marked region in water scenario dataset.

4.2.4 Experiments in Coastal Scenarios

Table 11 shows the good fusion performance of the proposed GP-based GO in coastal environments. GPbasedGO consistently outperforms or is highly competitive in all four test cases, demonstrating its strong ability to handle complex shoreline structures, mixed land-water boundaries, and significant spectral heterogeneity, all of which are common in coastal remote sensing scenes. These factors make coastal image fusion far more difficult than in homogeneous inland areas.

In QuickBird (Image 1), GPbasedGO outperforms all six evaluation metrics: ERGAS (6.092898), SAM (8.71684), RASE (26.92863), RMSE (25.08116), UIQI (0.954367), and CC (0.992877). Compared to the original GO, the suggested method greatly improves both spectral and structural consistency, particularly in UIQI and CC, indicating that GP-basedGO is more effective at retaining fine shoreline textures and decreasing spectral distortion in coastal transition zones.

In the tough QuickBird (Image 2) environment, GPbasedGO achieves the best values in RASE, RMSE, UIQI, and CC, while its ERGAS and SAM remain highly competitive. Although SA and DE perform marginally better in ERGAS and SAM, respectively, GP-basedGO performs much better overall in similarity-related metrics. In particular, as compared to GO, the suggested method increases UIQI from 0.928916 to 0.944103 and CC from 0.914426 to 0.938331, demonstrating its improved capacity to maintain structural continuity and edge integrity in extremely varied coastal regions.

GF (Image 1) and GF (Image 2) datasets show similar advantages. For GF (Image 1), GP-basedGO outperforms all six assessment parameters, demonstrating exceptional retention of fine coastline details and spectral fidelity. In GF (Image 2), the suggested technique achieves the best performance in all six metrics, with especially considerable gains in ERGAS, RASE, RMSE, UIQI, and CC when compared to competing algorithms. These findings demonstrate that the Gaussian-process-guided search mechanism can significantly improve optimization accuracy and strike a better balance between spatial detail augmentation and spectrum information preservation across a variety of sensor attributes and coastal distributions.

The visual comparisons in Fig. 12 support the quantitative results. The fusion images created by GPbasedGO show crisper shoreline contours, finer border transitions between land and water, and more natural tonal fluctuations in shallow-water and reflective areas. Compared to baseline techniques, the suggested methodology effectively reduces spectral distortion and prevents boundary blurring, which is typical in coastal fusion tasks. Furthermore, the retention of faint texture information in ports, tidal flats, and coastline vegetation regions has been improved. These visual observations align with the superior objective results given in Table 11, highlighting the robustness and usefulness of GP-basedGO for remote sensing image fusion in coastal environments.

Figure 12: Visual assessment for marked region in coast dataset.

5  Conclusions

In this paper, a Genetic Programming-based Metaheuristic framework (GP-MAs) is developed to address the constraints of standard metaheuristic algorithms imposed by hand-crafted updating strategies. The suggested system, which incorporates Genetic Programming (GP) into the search process, allows for the automatic evolution of symbolic update equations, boosting the adaptability and flexibility of metaheuristic optimization. Based on this architecture, a hybrid optimizer called GPbasedGO is created, with the Growth Optimizer (GO) as the baseline algorithm. Extensive studies on the CEC2022 benchmark suite show that GP-basedGO has competitive convergence behavior, optimization performance, and stability when compared to many representative state-of-the-art algorithms. Furthermore, the efficacy of GP-basedGO is assessed in multispectral and panchromatic image fusion tasks across urban, vegetation, water, and coastal situations, with the suggested approach outperforming many image quality measures such as ERGAS, SAM, RMSE, UIQI, and CC. The experimental results show that the proposed GP-MAs framework offers a flexible and partially automated approach to symbolic strategy evolution in metaheuristic optimization. This study validates GP-MAs with GO as a representative human-inspired optimizer, proving the viability of GP-assisted update rule evolution for representative optimization problems. However, the current study does not fully cover highly dynamic, limited, or non-stationary optimization situations, and further validation of swarm-based and physics-based optimizers remains a future task. Future research will expand the suggested framework to multi-objective, large-scale, and discrete optimization issues, as well as investigate its integration with deep learning and reinforcement learning approaches.

Although the proposed GP-MAs framework achieves competitive optimization performance and demonstrates its effectiveness in remote sensing image fusion, some limitations remain. The theoretical convergence analysis of the evolved symbolic update rules has not been fully established, and the computational cost of the offline GP evolution process requires further quantitative investigation. Future research will focus on developing more rigorous theoretical analysis methods, improving the efficiency of the GP evolution process, and extending the proposed framework to a wider range of metaheuristic algorithms and complex optimization scenarios.

Acknowledgement: None.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, Shu-Chuan Chu and Hongmei Yang; methodology, Jeng-Shyang Pan, Wenda Li, and Hongmei Yang; software, Wenda Li; validation, Jeng-Shyang Pan and Wenda Li; formal analysis, Wenda Li, Zhi-Gang Du, and Lingping Kong; investigation, Wenda Li; resources, Shu-Chuan Chu; data curation, Wenda Li and Shu-Chuan Chu; writing—original draft preparation, Wenda Li; writing—review and editing, Shu-Chuan Chu; visualization, Wenda Li; supervision, Jeng-Shyang Pan, Zhi-Gang Du, Hongmei Yang, and Lingping Kong; project administration, Jeng-Shyang Pan; funding acquisition, Not applicable. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are openly available in GitHlub at https://github.com/liwenmda/Metaheuristic-Algorithms.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Dokeroglu T, Sevinc E, Kucukyilmaz T, Cosar A. A survey on new generation metaheuristic algorithms. Comput Ind Eng. 2019;137(5):106040. doi:10.1016/j.cie.2019.106040. [Google Scholar] [CrossRef]

2. Rajwar K, Deep K, Das S. An exhaustive review of the metaheuristic algorithms for search and optimization: taxonomy, applications, and open challenges. Artif Intell Rev. 2023;56(11):13187–257. doi:10.1007/s10462-023-10470-y. [Google Scholar] [PubMed] [CrossRef]

3. Nguyen TT, Pan JS, Dao TK. An improved flower pollination algorithm for optimizing layouts of nodes in wireless sensor network. IEEE Access. 2019;7:75985–98. doi:10.1109/access.2019.2921721. [Google Scholar] [CrossRef]

4. Du ZG, Pan JS, Chu SC, Luo HJ, Hu P. Quasi-affine transformation evolutionary algorithm with communication schemes for application of RSSI in wireless sensor networks. IEEE Access. 2020;8:8583–94. doi:10.1109/access.2020.2964783. [Google Scholar] [CrossRef]

5. Zhou G, Xie Y, Lan H, Tian W, Buyya R, Wu K. Information interaction and partial growth-based multi-population growable genetic algorithm for multi-dimensional resources utilization optimization of cloud computing. Swarm Evol Comput. 2024;87(4):101575. doi:10.1016/j.swevo.2024.101575. [Google Scholar] [CrossRef]

6. Abasi AK, Aloqaily M, Guizani M, Ouni B. Metaheuristic algorithms for 6G wireless communications: recent advances and applications. Ad Hoc Netw. 2024;158(4):103474. doi:10.1016/j.adhoc.2024.103474. [Google Scholar] [CrossRef]

7. Wang L, Pan Z, Wang J. A review of reinforcement learning based intelligent optimization for manufacturing scheduling. Complex Syst Model Simul. 2021;1(4):257–70. doi:10.23919/csms.2021.0027. [Google Scholar] [CrossRef]

8. Talbi EG. Machine learning into metaheuristics: a survey and taxonomy. ACM Comput Surv. 2021;54(6):1–32. doi:10.1145/3459664. [Google Scholar] [CrossRef]

9. Wang X, Wu Y. An improved deep spiking neural network model based on bionic optimization algorithm. J Netw Intell. 2024;9(1):160–75. doi:10.1109/icaci.2012.6463299. [Google Scholar] [CrossRef]

10. Elsayed HA, Hameed MF, El Sherbiny MM. Image classification using decision tree classifier and features extraction using Hough transform and genetic algorithm. J Inf Hiding Multimed Signal Process. 2025;16(1):438–52. [Google Scholar]

11. Mitchell M. An introduction to genetic algorithms. Cambridge, MA, USA: MIT Press; 1998. [Google Scholar]

12. Pan JS, Chen J, Li J, Meng Z. Differential evolution with exponential crossover: a survey. Comput Sci Rev. 2026;62(4):100990. doi:10.1016/j.cosrev.2026.100990. [Google Scholar] [CrossRef]

13. Beyer HG, Schwefel HP. Evolution strategies—a comprehensive introduction. Nat Comput. 2002;1:3–52. [Google Scholar]

14. Simon D. Biogeography-based optimization. IEEE Trans Evol Comput. 2008;12(6):702–13. doi:10.1109/tevc.2008.919004. [Google Scholar] [CrossRef]

15. Shams B. Using network properties to evaluate targeted immunization algorithms. Netw Biol. 2014;4(3):74. [Google Scholar]

16. Kennedy J, Eberhart R. Particle swarm optimization. In: Proceedings of ICNN’95-International Conference on Neural Networks; 1995 Nov 27–Dec 1; Perth, WA, Australia. p. 1942–8. [Google Scholar]

17. Blum C. Ant colony optimization: introduction and recent trends. Phys Life Rev. 2005;2(4):353–73. [Google Scholar]

18. Karaboga D, Akay B. A comparative study of artificial bee colony algorithm. Appl Math Comput. 2009;214(1):108–32. doi:10.1016/j.amc.2009.03.090. [Google Scholar] [CrossRef]

19. Mirjalili S, Mirjalili SM, Lewis A. Grey wolf optimizer. Adv Eng Softw. 2014;69:46–61. doi:10.1016/j.advengsoft.2013.12.007. [Google Scholar] [CrossRef]

20. Mirjalili S, Lewis A. The whale optimization algorithm. Adv Eng Softw. 2016;95(12):51–67. doi:10.1016/j.advengsoft.2016.01.008. [Google Scholar] [CrossRef]

21. Abdollahzadeh B, Gharehchopogh FS, Khodadadi N, Mirjalili S. Mountain gazelle optimizer: a new nature-inspired metaheuristic algorithm for global optimization problems. Adv Eng Softw. 2022;174(3):103282. doi:10.1016/j.advengsoft.2022.103282. [Google Scholar] [CrossRef]

22. Bertsimas D, Tsitsiklis J. Simulated annealing. Stat Sci. 1993;8(1):10–5. doi:10.1214/ss/1177011077. [Google Scholar] [CrossRef]

23. Rashedi E, Nezamabadi-Pour H, Saryazdi S. GSA: a gravitational search algorithm. Inf Sci. 2009;179(13):2232–48. doi:10.1016/j.ins.2009.03.004. [Google Scholar] [CrossRef]

24. Mirjalili S. SCA: a sine cosine algorithm for solving optimization problems. Knowl-based Syst. 2016;96(63):120–33. doi:10.1016/j.knosys.2015.12.022. [Google Scholar] [CrossRef]

25. Mirjalili S, Mirjalili SM, Hatamlou A. Multi-verse optimizer: a nature-inspired algorithm for global optimization. Neural Comput Appl. 2016;27(2):495–513. doi:10.1007/s00521-015-1870-7. [Google Scholar] [CrossRef]

26. Rao RV, Savsani VJ, Vakharia DP. Teaching-learning-based optimization: a novel method for constrained mechanical design optimization problems. Comput-Aided Des. 2011;43(3):303–15. doi:10.1016/j.cad.2010.12.015. [Google Scholar] [CrossRef]

27. Geem ZW, Kim JH, Loganathan GV. A new heuristic optimization algorithm: harmony search. Simulation. 2001;76(2):60–8. doi:10.1177/003754970107600201. [Google Scholar] [CrossRef]

28. Satapathy S, Naik A. Social group optimization (SGOa new population evolutionary optimization technique. Complex & Intell Syst. 2016;2(3):173–203. doi:10.1007/s40747-016-0022-8. [Google Scholar] [CrossRef]

29. Zhang W, Pan K, Li S, Wang Y. Special forces algorithm: a novel meta-heuristic method for global optimization. Math Comput Simul. 2023;213:394–417. [Google Scholar]

30. Pan JS, Zhang LG, Wang RB, Snášel V, Chu SC. Gannet optimization algorithm: a new metaheuristic algorithm for solving engineering optimization problems. Math Comput Simul. 2022;202:343–73. [Google Scholar]

31. Wang J, Chen Y, Lu C, Heidari AA, Wu Z, Chen H. The status-based optimization: algorithm and comprehensive performance analysis. Neurocomputing. 2025;647:130603. [Google Scholar]

32. Liang Z, Wang Z, Mohamed AW. A novel hybrid algorithm based on improved marine predators algorithm and equilibrium optimizer for parameter extraction of solar photovoltaic models. Heliyon. 2024;10(19):e38412. doi:10.1016/j.heliyon.2024.e38412. [Google Scholar] [PubMed] [CrossRef]

33. Song H, Bei J, Zhang H, Wang J, Zhang P. Hybrid algorithm of differential evolution and flower pollination for global optimization problems. Expert Syst Appl. 2024;237(14):121402. doi:10.1016/j.eswa.2023.121402. [Google Scholar] [CrossRef]

34. Zhao J, Chen D, Xiao R, Chen J, Pan JS, Cui ZH, et al. Multi-objective firefly algorithm with adaptive region division. Appl Soft Comput. 2023;147(7):110796. doi:10.1016/j.asoc.2023.110796. [Google Scholar] [CrossRef]

35. Ma Z, Guo H, Chen J, Li Z, Peng G, Gong YJ, et al. Metabox: a benchmark platform for meta-black-box optimization with reinforcement learning. Adv Neural Inf Process Syst. 2023;36:10775–95. doi:10.5555/3666122.3666596. [Google Scholar] [CrossRef]

36. Finn C, Abbeel P, Levine S. Model-agnostic meta-learning for fast adaptation of deep networks. In: Proceedings of the 34th International Conference on Machine Learning; 2017 Aug 6–11; Sydney, NSW, Australia. p. 1126–35. [Google Scholar]

37. Liu F, Tong X, Yuan M, Lin X, Luo F, Wang Z, et al. Evolution of heuristics: towards efficient automatic algorithm design using large language model. arXiv:2401.02051. 2024. [Google Scholar]

38. Drugan MM. Reinforcement learning versus evolutionary computation: a survey on hybrid algorithms. Swarm Evol Comput. 2019;44:228–46. [Google Scholar]

39. Salimans T, Ho J, Chen X, Sidor S, Sutskever I. Evolution strategies as a scalable alternative to reinforcement learning. arXiv:1703.03864. 2017. [Google Scholar]

40. Guo H, Ma Y, Ma Z, Chen J, Zhang X, Cao Z, et al. Deep reinforcement learning for dynamic algorithm selection: a proof-of-principle study on differential evolution. IEEE Trans Syst Man Cybern Syst. 2024;54(7):4247–59. [Google Scholar]

41. Ma Z, Guo H, Chen J, Peng G, Cao Z, Ma Y, et al. LLaMoCo: instruction tuning of large language models for optimization code generation. arXiv:2403.01131. 2024. [Google Scholar]

42. Gomes HS, Léger B, Gagné C. Meta learning black-box population-based optimizers. arXiv:2103.03526. 2021. [Google Scholar]

43. AhmadiTeshnizi A, Gao W, Udell M. OptiMUS: scalable optimization modeling with (MI) LP solvers and large language models. arXiv:2402.10172. 2024. [Google Scholar]

44. Wang J, Zheng Y, Zhang Z, Peng H, Wang H. A novel multi-state reinforcement learning-based multi-objective evolutionary algorithm. Inf Sci. 2025;688(9):121397. doi:10.1016/j.ins.2024.121397. [Google Scholar] [CrossRef]

45. Ma Z, Guo H, Gong YJ, Zhang J, Tan KC. Toward automated algorithm design: a survey and practical guide to meta-black-box-optimization. IEEE Trans Evol Comput. 2025;30(2):667–87. [Google Scholar]

46. Guo H, Ma S, Huang Z, Hu Y, Ma Z, Zhang X, et al. Reinforcement learning-based self-adaptive differential evolution through automated landscape feature learning. In: Proceedings of the Genetic and Evolutionary Computation Conference; 2025 Jul 14–18; Malaga, Spain. p. 1117–26. [Google Scholar]

47. Koza JR. Genetic programming as a means for programming computers by natural selection. Stat Comput. 1994;4(2):87–112. doi:10.1007/bf00175355. [Google Scholar] [CrossRef]

48. Xue X, Lin JCW. Heterogeneous sensor integration through co-evolutionary genetic programming with multiple individual representations. Swarm Evol Comput. 2025;98(5):102098. doi:10.1016/j.swevo.2025.102098. [Google Scholar] [CrossRef]

49. Zhang Q, Gao H, Zhan ZH, Li J, Zhang H. Growth optimizer: a powerful metaheuristic algorithm for solving continuous and discrete global optimization problems. Knowl Based Syst. 2023;261:110206. [Google Scholar]

50. Peng Y, Gu S, Liang Y, Ouyang K, Li Y, Wang K, et al. Wave optics optimizer: a novel meta-heuristic algorithm for engineering optimization. Commun Nonlinear Sci Numer Simul. 2026;152:109337. [Google Scholar]

51. Abdullah M, Alimgeer KS, Hafeez G, Alghamdi B, Alsafran AS, Babar MZ. Electromagnetic wave propagation algorithm: a novel electromagnetic wave propagation-inspired optimizer for engineering applications. Ain Shams Eng J. 2025;16(11):103615. [Google Scholar]

52. Taheri A, RahimiZadeh K, Baumbach J, Beheshti A, Zolotareva O, Al-Betar MA, et al. Farthest better or nearest worse optimizer: a novel metaheuristic algorithm. Artif Intell Rev. 2025;59(2):79. doi:10.1007/s10462-025-11443-z. [Google Scholar] [CrossRef]

53. Sharifi T, Mirsalim M, Soleimanian Gharehchopogh F, Mirjalili S. Cultural history optimization algorithm: a new human-inspired metaheuristic algorithm for engineering optimization problems. Neural Comput Appl. 2025;37(25):21009–68. doi:10.1007/s00521-025-11379-z. [Google Scholar] [CrossRef]

54. Chen H, Ahmadianfar I, Heidari AA, Kordani M, Samadi Koucheksaraee A, Liang G. LEE: a physics-inspired optimizer based on LangEvin equation. Neurocomputing. 2026;666:132288. [Google Scholar]