Paper deep dive
SeisEvo: Evolution of Seismic Data Reconstruction Algorithms by Agents
Yingjie Xu, Siwei Yu, Jianwei Ma
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 8/20/2026, 4:14:07 AM
Summary
The paper introduces SeisEvo, a framework for the autonomous discovery of seismic data reconstruction algorithms using LLM-driven multi-agent program evolution. Unlike deep learning methods that rely on opaque neural weights, SeisEvo searches for explicit, white-box algorithms by modifying classical seed algorithms (POCS and MSSA) under strict physical and mathematical constraints. The system discovered Evo-POCS, which improves SNR by 3.49 dB over classic POCS, and Evo-MSSA, which improves SNR by >7 dB over classic MSSA and >3 dB over ODRR, demonstrating that agentic evolution can complement deep learning by producing inspectable and deployable operators.
Entities (9)
Relation Signals (7)
SeisEvo → discovers → Evo-MSSA
confidence 95% · For simultaneous interpolation and denoising, it discovered a reliability-grouped singular-value shrinkage; Evo-MSSA improves the average reconstruction SNR by more than 7 dB over classic MSSA
SeisEvo → discovers → Evo-POCS
confidence 95% · For interpolation without added noise, the search discovered a residual-gated, phase-aligned dip-consistency projection; Evo-POCS improves the SNR over classic POCS by 3.49 dB
Evo-POCS → improves → POCS
confidence 95% · Evo-POCS improves the SNR over classic POCS by 3.49 dB on average across missing ratios from 30% to 70%.
Evo-MSSA → improves → MSSA
confidence 95% · Evo-MSSA improves the average reconstruction SNR by more than 7 dB over classic MSSA
SeisEvo → isbasedon → LLM
confidence 95% · an LLM-driven multi-agent search modifies only the components that the user has opened for editing
SeisEvo → uses → EVE
confidence 95% · SeisEvo uses the Evolutionary Ensemble of Agents (EvE) (38), a decentralized evolutionary framework that organizes multiple coding agents to iteratively propose and refine candidate programs.
Evo-MSSA → improves → ODRR
confidence 90% · Evo-MSSA improves the average reconstruction SNR by more than 7 dB over classic MSSA and by more than 3 dB over a stronger rank-reduction baseline.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Classical seismic data reconstruction relies on manually designed structural priors and iterative operators, whose coupled design space is far larger than manual trial and error can explore systematically. Deep-learning methods encode the reconstruction rules in learned weights rather than in an explicit operator that can be inspected and modified. We propose SeisEvo (Seismic Algorithm Evolution), which does not optimize a single reconstruction result but searches for the algorithm that produces it. Starting from a classical reconstruction algorithm, an LLM-driven multi-agent search modifies only the components that the user has opened for editing, without prescribing the mechanism to be discovered. Candidates that violate the physical constraints of the task are rejected outright, and the remaining ones are scored by execution. The output is neither an agent system nor a neural network, but a standalone white-box algorithm that requires no agent or neural network at inference time. For interpolation without added noise, the search discovered a residual-gated, phase-aligned dip-consistency projection; Evo-POCS improves the SNR over classic POCS by 3.49 dB on average across missing ratios from 30% to 70%. For simultaneous interpolation and denoising, it discovered a reliability-grouped singular-value shrinkage; Evo-MSSA improves the average reconstruction SNR by more than 7 dB over classic MSSA and by more than 3 dB over a stronger rank-reduction baseline. Both operators retain their gains on data not used during the search. To the best of our knowledge, this is the first study to formulate the design of a seismic reconstruction operator as a constrained, LLM-driven program evolution task. Agentic algorithm evolution can thus complement deep learning in discovering explicit, inspectable, and deployable seismic processing algorithms.
Tags
Links
- Source: https://arxiv.org/abs/2608.18272v1
- Canonical: https://arxiv.org/abs/2608.18272v1
Trouble viewing inline? Open PDF directly →
Full Text
92,385 characters extracted from source content.
Expand or collapse full text
SeisEvo: Evolution of Seismic Data Reconstruction Algorithms by Agents Yingjie Xu Affiliation: School of Mathematics and Center of Geophysics, Harbin Institute of Technology, Harbin, Heilongjiang, China Affiliation: Department of Mathematics, National University of Singapore, Singapore Siwei Yu Thanks: Corresponding author: siweiyu@hit.edu.cn Affiliation: School of Mathematics and Center of Geophysics, Harbin Institute of Technology, Harbin, Heilongjiang, China Jianwei Ma Affiliation: School of Mathematics and Center of Geophysics, Harbin Institute of Technology, Harbin, Heilongjiang, China Affiliation: School of Earth and Space Sciences, Peking University, Beijing, China Abstract Classical seismic data reconstruction relies on manually designed structural priors and iterative operators, whose coupled design space is far larger than manual trial and error can explore systematically. Deep-learning methods encode the reconstruction rules in learned weights rather than in an explicit operator that can be inspected and modified. We propose SeisEvo (Seismic Algorithm Evolution), which does not optimize a single reconstruction result but searches for the algorithm that produces it. Starting from a classical reconstruction algorithm, an LLM-driven multi-agent search modifies only the components that the user has opened for editing, without prescribing the mechanism to be discovered. Candidates that violate the physical constraints of the task are rejected outright, and the remaining ones are scored by execution. The output is neither an agent system nor a neural network, but a standalone white-box algorithm that requires no agent or neural network at inference time. For interpolation without added noise, the search discovered a residual-gated, phase-aligned dip-consistency projection; Evo-POCS improves the SNR over classic POCS by 3.49 dB on average across missing ratios from 30% to 70%. For simultaneous interpolation and denoising, it discovered a reliability-grouped singular-value shrinkage; Evo-MSSA improves the average reconstruction SNR by more than 7 dB over classic MSSA and by more than 3 dB over a stronger rank-reduction baseline. Both operators retain their gains on data not used during the search. To the best of our knowledge, this is the first study to formulate the design of a seismic reconstruction operator as a constrained, LLM-driven program evolution task. Agentic algorithm evolution can thus complement deep learning in discovering explicit, inspectable, and deployable seismic processing algorithms. Keywords: Seismic data reconstruction; Algorithm evolution; Large language model agents; Interpretable operators 1 Introduction Seismic data reconstruction is an important processing step for improving the spatial sampling quality of seismic records and for restoring the continuity of subsurface reflection events. Because of surface obstacles, economic cost, and practical acquisition conditions, field seismic data commonly suffer from missing traces and irregular spatial sampling, which degrade the reliability of subsequent imaging, inversion, and interpretation (13; 28). Since the missing samples are not observed, seismic data reconstruction is an underdetermined inverse problem, and prior information about seismic structure must be introduced to constrain the solution, such as transform-domain sparsity (2; 13), local coherence (26; 21), and structured low-rank behavior (23). Although data-driven methods have advanced rapidly in recent years, many seismic reconstruction methods, as in many other seismic processing steps, still rely on explicitly designed mathematical operators and structural priors: how a prior is represented, through which computational operator it is imposed, and how the estimate is kept consistent with the observations are usually decided manually from domain experience. Reconstruction performance therefore depends largely on how effectively these operators encode seismic structure. Classical reconstruction methods encode these structural priors explicitly as interpretable algorithmic procedures. Based on transform-domain sparsity, the projection onto convex sets (POCS) method (2) recovers missing traces by alternating transform-domain thresholding with data-consistency enforcement. Based on local coherence, prediction-filter and local-slope-based methods (26; 21; 33) instead extrapolate missing traces directly from the predictable relation between neighboring traces. Based on structured low-rank behavior, multichannel singular spectrum analysis (MSSA) (23) applies a truncated low-rank projection to block-Hankel embeddings, restoring the low-rank structure formed by coherent reflection events and suppressing the rank increase caused by missing traces and noise. These methods rest on explicit structural assumptions, follow transparent computational steps, and are easy to deploy. POCS with exponential hard thresholding (10) and block-Hankel MSSA with truncated singular value decomposition (23) are widely used as reference implementations of the two families, and are adopted as the seed algorithms in this work. Building on these classical templates, many manually designed improvements have been proposed. For POCS, they mainly concern the choice of the sparsifying transform (29; 40), the form of the thresholding or shrinkage function, and the threshold decay schedule used during the iterations (10; 11). For MSSA, they mainly concern the embedding and local window design, the rank selection, and the singular-value processing rule (39; 31), leading further to enhanced low-rank methods such as damped rank reduction (14; 6) and optimally damped rank reduction (ODRR) (5). These improvements are introduced through expert-guided, mechanism-specific design and are usually validated under specific acquisition conditions or noise levels. However, most studies address only one or a few dimensions of the design space, and their effectiveness depends on the particular data distribution, sampling density, and signal-to-noise ratio, so they may require retuning or redesign under different conditions. Moreover, these design factors interact: the threshold form and the decay schedule, or the rank selection and the singular-value shrinkage, jointly determine performance, so adjusting each factor in isolation is not sufficient for exploring their combined effects. The design space left open by these classical templates is therefore far larger than manual exploration can systematically cover. In recent years, deep learning has provided an alternative route for seismic data reconstruction. By learning data-driven reconstruction mappings or priors, such methods have been applied to seismic data interpolation and noise suppression (36; 30; 15; 7). To reduce the reliance on paired complete data, unsupervised and self-supervised methods further construct training signals from unlabeled observations or from the data to be reconstructed themselves (17; 20; 1; 4; 32). Seismic foundation models have also emerged, which learn general seismic representations through large-scale self-supervised pretraining and transfer them to different downstream tasks such as interpolation, denoising, and inversion (25; 8). These methods differ in nature from classical reconstruction algorithms. The performance of supervised methods depends on the data distribution covered by the training set, and their generalization needs to be re-examined when the acquisition conditions, noise levels, or geological structures differ appreciably from the training data. Unsupervised and self-supervised methods reduce the reliance on fully labeled data, but their performance is still affected by the network architecture, the regularization, and the training strategy. For most end-to-end methods, the main reconstruction rules are distributed across the network weights and the training procedure, and are therefore difficult to write down, inspect, and modify under given constraints in the way classical operators are. Deep learning has thus not removed manual design; it has moved part of the design effort from the reconstruction operator to the network architecture, the loss function, and the training strategy. Classical and deep-learning methods still share a common feature: both optimize within a computational structure that is defined in advance by humans. In most existing studies on seismic reconstruction, the task-relevant algorithmic structure is still designed manually, and treating the algorithmic structure itself as the object of search remains rare. In recent years, large language model (LLM)-driven program search and evolution have opened a new possibility for algorithm design. These methods represent a candidate algorithm as executable code, let an LLM propose modifications, and then run and score each candidate with an automated evaluator, using the resulting feedback to guide the subsequent search. FunSearch showed that this paradigm can discover improved mathematical constructions and heuristics in combinatorial mathematics and combinatorial optimization (24). EoH and ReEvo subsequently improved heuristic generation, the use of feedback, and search efficiency, although the object of evolution there is typically a heuristic rule inside a fixed algorithmic pipeline (16; 35). In an applied setting, Eureka evolved reward code and obtained reward functions that outperform human-engineered ones on most robotic control tasks (18). AlphaEvolve further extended the paradigm to entire code bases and more general algorithm optimization tasks, and applied it to problems such as matrix multiplication, data center scheduling, and hardware design (22). Unlike end-to-end neural approaches, such a search returns a program that can be executed and inspected on its own rather than a set of learned network weights; interpretability, however, is not obtained automatically, but depends on how the program is represented and on the constraints imposed during the search. These advances indicate that the algorithm itself can become the object of search. Seismic data reconstruction offers a well-suited setting for this paradigm. First, reconstruction quality can be measured directly by metrics such as SNR and SSIM on the prescribed search-time blocks for which complete references are available, so that each candidate can be executed and assigned an explicit score. Second, classical reconstruction algorithms have a clear and modular structure, so that mechanisms such as the thresholding rule, the decay schedule, and the low-rank projection can be designated individually as editable components. Third, seismic reconstruction carries admissibility conditions that are explicit and decidable: how the observations are used, whether the sampling mask is left unchanged, and whether the output is real-valued are all determined once the task is given, and can therefore be checked programmatically. Applying this paradigm to seismic reconstruction, however, is not a matter of simply running a generic search. When domain legality conditions are absent, an increase in score does not necessarily correspond to an improvement of the algorithm: a candidate may reach a physically meaningless score by corrupting the observations or altering the sampling mask; it may raise a metric through unexplained free constants or dataset-specific shortcuts rather than by finding a new mechanism. For seismic reconstruction, legality depends not only on mathematical and numerical correctness but also on physical measurement conditions, so how these conditions are formalized and enforced as part of the search is itself part of the problem. To date, the use of LLM-driven program evolution to systematically discover standalone, white-box seismic reconstruction operators remains largely unexplored, and validation across different mathematical mechanisms and different reconstruction task settings is especially lacking. To address these issues, we propose SeisEvo (Seismic Algorithm Evolution), a domain-grounded framework for the autonomous discovery of seismic reconstruction algorithms. SeisEvo takes the reconstruction algorithm itself, rather than a particular reconstructed result, as the object of search: instead of optimizing the output of a single reconstruction, it searches for an algorithm that produces high-quality reconstructions. A run is specified by a domain expert through a transparent classical seed algorithm A0A_0, an editable operator surface S, a reconstruction-quality objective O, a constraint set C, and a search budget T. Here S states which mechanisms may be modified, such as the thresholding rule, the decay schedule, and the low-rank projection, whereas the task-specific data-fidelity requirements remain constrained throughout the search. The expert thus does not specify the final form of the operator in advance, but defines the search boundary, the physical contracts, and the scoring rules, leaving the mechanism itself to be found by the search. SeisEvo uses the Evolutionary Ensemble of Agents (EvE) (38), a decentralized evolutionary framework that organizes multiple coding agents to iteratively propose and refine candidate programs. EvE and EvE-based workflows have been used for scientific problems including fluid-control design, operator learning, and PDE symbolic discovery (27; 34; 37). Here it serves as the generic multi-agent search backend. Our contribution lies in the domain-specific framework for seismic algorithm evolution, including the algorithm representation, the editable operator surface, the legality constraints, the execution-based scoring, and the external validation protocol. When the search ends, the highest-scoring legal candidate is returned as A⋆A , a standalone white-box reconstruction operator. Its computational steps and introduced parameters can be written down, inspected, and reproduced, and it requires no agent, prompt, or neural network at inference time. We instantiate SeisEvo on two classical algorithm families with different mathematical foundations: POCS based on Fourier-domain sparsity and MSSA based on block-Hankel low-rank structure. The resulting operators are referred to as Evo-POCS and Evo-MSSA, respectively. This paper makes three main contributions. First, we propose SeisEvo, an algorithm evolution framework for seismic data reconstruction that formulates the design of a reconstruction algorithm as a domain-constrained autonomous discovery task. Starting from a classical seed algorithm, SeisEvo searches for new operator mechanisms that improve reconstruction quality, while task-specific data-fidelity requirements remain enforced throughout the search. Second, we establish an executable search and validation protocol that uses hard legality gates and requires mechanism-level changes with traceable parameter origins. External data are accessed only after the search, once the selected operator has been fixed. Third, we validate SeisEvo on two classical algorithm families with different mathematical foundations. For interpolation of a field dataset without added synthetic noise, Evo-POCS improves the SNR over classic POCS by 3.49 dB on average across missing ratios from 30% to 70%. For simultaneous interpolation and denoising, Evo-MSSA improves the mean reconstruction SNR by more than 7 dB over classic MSSA and by more than 3 dB over the ODRR baseline over the tested noise range. Both evolved operators retain their performance gains on synthetic and field data that were not used during the search. To the best of our knowledge, this is the first reported application of LLM-driven constrained program evolution to the discovery of standalone seismic reconstruction operators. The remainder of this paper is organized as follows. Section 2 introduces the seismic reconstruction problem and the SeisEvo framework, including the task specification and constraint design. Section 3 presents the POCS and MSSA case studies and their experimental results. Section 4 discusses what the search contributed, the scope and cost of the protocol, and the open problems it leaves. Section 5 concludes the paper. 2 Method 2.1 Seismic data reconstruction problem Let X denote a complete seismic volume of size nt×nx×nyn_t× n_x× n_y, and let M be a binary sampling mask of the same size, where M(t,x,y)=1M(t,x,y)=1 indicates an acquired sample and M(t,x,y)=0M(t,x,y)=0 a missing one. Since traces are missing as whole columns, M is constant along the time axis. Let N denote additive noise on the acquired samples. The observed incomplete volume is Y=M⊙(X+N),Y=M (X+N), (1) where ⊙ denotes element-wise multiplication, and N=0N=0 gives the noise-free interpolation setting. The task is to estimate a complete volume X X from (Y,M)(Y,M). Because the missing traces are unobserved, the problem is underdetermined, and it is commonly formulated as a regularized inverse problem, X^=argminZ12‖Y−M⊙Z‖F2+λR(Z), X= _Z\; 12\, \|Y-M Z \|_F^2\;+\;λ\,R(Z), (2) where the first term enforces fidelity to the acquired traces and R(Z)R(Z) is a regularization term encoding prior knowledge about the structure of seismic data, balanced by λ. Two representative and widely used priors are transform-domain sparsity and structured low rank: the former assumes that seismic data admit compact representations in a suitable transform domain, while the latter assumes that coherent seismic events give rise to approximately low-rank structures after suitable embedding, which missing traces and noise tend to break. Problem (2) motivates a broad class of iterative schemes that alternate between enforcing the assumed prior and honoring the observations, Z(j+1)=data(prior(j)(Z(j))),Z^(j+1)=P_data\! (T_prior^(j)(Z^(j)) ), (3) where prior(j)T_prior^(j) is an iteration-dependent operator enforcing the seismic prior, and dataP_data is the task-specified data-fidelity step. In noise-free interpolation it is the hard projection data(U)=Y+(1−M)⊙U,P_data(U)=Y+(1-M) U, (4) which restores the acquired traces exactly; when the acquired samples are noisy, the hard projection is replaced by a task-specific relaxed fidelity step that balances the observations and the current prediction. Its concrete form is given in the corresponding case study. POCS and MSSA are two representative classical solvers associated with these priors and serve as the seed algorithms in this work: POCS realizes prior(j)T_prior^(j) through transform-domain thresholding, encoding the sparsity prior, while MSSA realizes it through a truncated low-rank projection of block-Hankel embeddings, encoding the low-rank prior. Taking these two algorithms as templates, numerous studies have developed further improvements. Improvements to POCS mostly focus on the sparsifying transform, the thresholding or shrinkage function, and the decay schedule (10; 29; 40; 9), with hard thresholding and exponential decay as the standard baseline. Improvements to MSSA mostly focus on the rank-reduction operator and rank selection (14; 39; 3), with truncated SVD and a fixed rank as the standard baseline. Such designs are interpretable and often effective, but each modification explores only a small portion of the possible design choices: tuning relies heavily on expert experience, and a setting that works under one sampling or noise condition rarely transfers unchanged to another. The design space left open by these templates is far larger than what manual exploration can cover. The recent emergence of LLM-based coding agents makes it feasible to explore this design space automatically. Motivated by this, we propose SeisEvo, which turns the open design space into a search space while keeping the mask-based structure of the task-specific data-fidelity step and the physical-legality constraints fixed. We instantiate it on the two seed algorithms introduced above: POCS with the sparsity prior and MSSA with the low-rank prior. These two solvers are mechanistically distinct, which lets us verify that the framework operates beyond a single algorithm family. The framework is detailed next. 2.2 The SeisEvo framework We introduce SeisEvo (Seismic Algorithm Evolution), a framework for the autonomous discovery of interpretable seismic reconstruction algorithms. SeisEvo treats the reconstruction algorithm itself, rather than a particular reconstructed result, as the object of search. Given a transparent classical seed algorithm A0A_0 of the form (3), a SeisEvo run with round budget T is written as A⋆=SeisEvo(A0,,O,,T),A =SeisEvo (A_0,\;S,\;O,\;C;\;T ), (5) where S denotes the editable parts of the algorithm, such as the thresholding function, the decay schedule, or the low-rank projection, while the mask-based structure of the task-specific data-fidelity step remains fixed. O is the objective, a measured reconstruction-quality score, and C is the constraint set, covering physical-legality constraints, parameter-provenance rules, anti-leakage requirements, and the search-time scoring protocol. The budget T is the number of search rounds. The returned A⋆A is the highest-scoring legal candidate, namely the evolved white-box algorithm. The defining principle is that the task specification states what counts as a legal and worthwhile algorithm, not what the algorithm should be. The objective O and constraint set C determine only how candidates are admitted and scored. They contain no candidate code, no target operator, and no fixed implementation template, so the mechanism itself is left for the search to discover. The output is a standalone white-box algorithm whose computational steps and introduced parameters can be inspected and executed without agents, prompts, or neural networks at inference time. A SeisEvo run consists of four stages: task instantiation, candidate generation, gating and scoring, and evidence binding and context update (Fig. 1). Stage 1 is performed once by a human. Stages 2 to 4 form a closed loop that repeats for T rounds. The four stages answer four questions in order: what algorithm is worth searching for, what to try next, whether a candidate qualifies to be compared, and what a completed trial counts as evidence for. External validation is not part of the loop; it is performed once on the final algorithm after the search ends. Stage 1: Task instantiation (human). A human defines what counts as a legal and worthwhile algorithm, not what the algorithm should be. This includes the seed algorithm A0A_0, the editable components S, the objective O, and the constraint set C defined above. Stage 2: Candidate generation (agents). Multiple coding agents work in parallel. Each agent reads sampled high-scoring candidates, failed candidates, diverse candidates, and the current search context. It then proposes new candidates by editing the algorithm only within S. Agents may also propose updates to the search context for later rounds, but these are treated as prospective hypotheses rather than adjudicated evidence. To support attribution and inspection, every candidate carries a mandatory rationale: what mechanism it changes, where each new parameter comes from, and why the change may help. Stage 3: Gating and scoring (framework). Physical legality is a hard gate, not a soft penalty. Candidates are rejected outright if they violate the task-specified data-consistency requirements, for example by breaking observed traces or altering the sampling mask in interpolation, or if they fail to run or diverge numerically. Surviving candidates are run on the search-time sub-blocks, scored by a common metric, and then admitted to the candidate pool. Rejected and failed candidates are kept as negative examples. Stage 4: Evidence binding and context update (framework). The framework records each candidate’s rationale, score, or failure log. It then updates which candidates and mechanism rationales are more likely to be shown to agents in the next round. This stage produces no new algorithmic hypotheses and contains no language-model judgment; it applies fixed programmatic rules to execution-derived scores and failure records. Agents only propose; all adjudication is done by the framework from real execution. After T rounds, the highest-scoring legal candidate under O is frozen as the evolved algorithm A⋆A , which then proceeds to external validation outside the search. Figure 1: The SeisEvo framework. A one-time task instantiation (seed algorithm and specification) is followed by a closed loop of candidate generation, gating and scoring, and evidence binding with context update, repeated for T rounds. The best candidate is returned as a standalone white-box operator. External validation, including full-volume generalization and robustness tests under different sampling and noise conditions, is performed once after the search and is not part of the loop. Algorithm 1 SeisEvo design run 1: seed algorithm A0A_0; editable surface S; objective O; constraint set C; rounds T 2: evolved white-box algorithm A⋆A 3: Stage 1 (human): instantiate (A0,,O,)(A_0,S,O,C) 4: pool ←A0←\A_0\ 5: for t=1,…,Tt=1,…,T do 6: Stage 2 (agents): agents read sampled candidates, failure records, and context, then propose new candidates by editing only within S, each with a rationale and parameter provenance 7: Stage 3 (framework): reject candidates that violate C, fail to run, or diverge; otherwise run them on search-time sub-blocks, score them by O, and add them to the pool 8: Stage 4 (framework): record each candidate’s rationale and outcome; update next-round context and read priorities by fixed programmatic rules 9: end for 10: return the best-scoring candidate A⋆A under O 2.3 Task specification Besides the seed algorithm A0A_0 and the editable surface S, a SeisEvo run is specified by the objective O and the constraint set C. Together they state what counts as a valid and worthwhile algorithm, without prescribing the mechanism to be found. Objective O. A scalar reconstruction-quality score. It may be a single metric, such as the signal-to-noise ratio (SNR), or a weighted combination of several metrics, such as SNR together with SSIM or MSE. A composite scalar objective may guide the search toward more robust results than a single metric alone (22). The concrete scoring design is task-dependent and is given in the experiments. Physical legality. For the interpolation tasks considered here, observed traces must be preserved exactly; for noisy reconstruction tasks, this hard projection can be replaced by a noise-aware fidelity rule. In addition, every candidate must leave the sampling mask unchanged, return a real-valued volume, and keep the fixed input–output interface. These are hard gates: violating any one rejects the candidate outright. Dataset-specific shortcuts and any use of test data during the search are likewise forbidden. Interpretability rules. Two rules keep the discovered algorithm auditable. First, no pure parameter tuning: the specification prohibits candidates whose only claimed change is retuning existing scalar constants, for example making a threshold decay faster. Changing the functional form of a threshold, schedule, or projection counts as a mechanism-level modification. Second, parameter provenance: every new free parameter must be inherited from the seed protocol, derived from observed statistics or mathematical bounds, introduced solely for numerical stability, or selected from a predeclared finite structural set. Unexplained free constants are forbidden. The first rule prevents progress by tuning numbers; the second requires an origin for any number that does change. Scoring during the search. Candidates are scored only on search-time sub-blocks extracted from the data; the external test volumes are never touched until the search has ended. These rules form the common basis shared by both case studies. The framework does not fix the exact form of the constraints: for a different task, the expert may add further physical or numerical constraints, such as preserving a given frequency band or bounding amplitude changes, which act as hard gates in the same way. When the search ends, the highest-scoring legal candidate A⋆A is frozen as the evolved algorithm: a standalone white-box solver that runs without any agent, prompt, or neural network at inference time. Because candidates are scored only on search-time sub-blocks during the search, we validate A⋆A once after the search to examine whether its performance generalizes beyond the data used to guide the search. In the experiments, this validation combines full-volume generalization and robustness tests under different sampling and noise conditions. We instantiate the protocol on POCS and MSSA, yielding Evo-POCS and Evo-MSSA. 3 Numerical experiments 3.1 Experimental protocol The data used for evolution are 64×64×6464× 64× 64 sub-blocks extracted from complete 3D seismic volumes. During the search, candidate algorithms are generated, scored, and used to update the context only on these sub-blocks. After the search ends, the discovered algorithm is tested on seismic data outside this set to assess its generalization. Figure 2: Search-time seismic sub-blocks extracted from complete 3D volumes. In seismic data reconstruction, the signal-to-noise ratio (SNR) is a common measure of reconstruction quality; we report SNR together with the structural similarity (SSIM) and the mean squared error (MSE): SNR=10log10‖X‖F2‖X−X^‖F2,SNR=10 _10 \|X\|_F^2\|X- X\|_F^2, (6) MSE=1Nsam‖X−X^‖F2,MSE= 1N_sam\|X- X\|_F^2, (7) SSIM(X,X^)=(2μXμX^+c1)(2σXX^+c2)(μX2+μX^2+c1)(σX2+σX^2+c2),SSIM(X, X)= (2 _X _ X+c_1)(2 _X X+c_2)( _X^2+ _ X^2+c_1)( _X^2+ _ X^2+c_2), (8) where X is the original complete data and X X the reconstruction, NsamN_sam is the total number of samples, and μ and σ denote means and (co)variances, and c1,c2c_1,c_2 are small constants that avoid a zero denominator. A higher SNR, a lower MSE, and an SSIM closer to one indicate better reconstruction. For multiple volumes or cases, each metric is computed per case and then averaged. 3.2 SeisEvo test on POCS for interpolation Task instantiation. The first case study instantiates SeisEvo on the POCS family for interpolation of randomly missing traces without adding synthetic noise. In this setting, no additional noise term is introduced in Eq. (1), i.e., N=0N=0, so the observation reduces to Y=M⊙XY=M X, where M is a random spatial sampling mask that drops whole traces. The four inputs of the method are instantiated as follows. The seed A0A_0 is classic exponential hard-threshold POCS. The editable surface S contains the per-iteration thresholding and shrinkage rules, the decay schedule, and any added projection or correction operator, but not the input–output interface or the observed-trace restoration. The constraint set C requires the observed traces to be restored exactly, the sampling mask to be preserved, and the output to be real-valued. It also forbids test-set leakage and unexplained free parameters. The objective is OPOCS=0.7clip(SNR18,0,1)+0.3SSIM,O_POCS=0.7\,clip\! ( SNR18,0,1 )+0.3\,SSIM, (9) where 18 dB is an empirical reference used during the search to scale the SNR term to a range comparable with SSIM. Classic POCS seed. The seed is classic POCS (2), with the exponential threshold decay widely adopted in later work (10). After a temporal Fourier transform, each frequency slice is reconstructed by alternating a sparse projection in the spatial-frequency domain with an exact observed-trace projection in the data domain. For a frequency slice Uk(j)U_k^(j), one iteration is Vk(j)=ℱ−1[Hτk,j(ℱUk(j))],V_k^(j)=F^-1 [H_ _k,j (F\,U_k^(j) ) ], (10) Uk(j+1)=Yk+(1−M)⊙Vk(j),U_k^(j+1)=Y_k+(1-M) V_k^(j), (11) where ℱF is the 2D spatial Fourier transform over a frequency slice, HτH_τ is hard thresholding, and YkY_k is the observed frequency slice. The threshold τk,j _k,j decays exponentially from 90% to 1% of the maximum transform-domain amplitude over 30 reconstruction iterations. The same threshold schedule and iteration count are used by classic POCS and Evo-POCS. Equation (11) is the data-consistency projection of Eq. (4) applied per frequency slice. This seed exploits Fourier-domain sparsity but does not explicitly use other properties of seismic data. Its reconstruction is therefore limited when dipping or laterally continuous events are not well represented by Fourier thresholding alone. The next part describes the addition that SeisEvo discovers on top of this seed. Search process. Starting from the classic POCS seed, SeisEvo runs for 20 rounds. Figure 3 shows the search-time objective OPOCSO_POCS as a function of the search round. The objective rises quickly from 0.713 at the seed within the first ten rounds, reaches its best value of 0.909 at round 14, and stays unchanged over the remaining six rounds, indicating that the 20-round budget is sufficient for this task. Correspondingly, the SNR on the training sub-blocks improves from 11.32 to 15.88 dB and SSIM from 0.911 to 0.972. This gain comes from a new structural step introduced by the SeisEvo search, not from retuning the threshold or decay parameters. Figure 3: Search-time objective OPOCSO_POCS over 20 search rounds. The score rises from 0.713 at the classic POCS seed to 0.909 at round 14 and remains unchanged afterwards. Evo-POCS. The new structural step introduced by the search is a data-driven local dip-consistency correction, applied after each Fourier thresholding step. In one sentence, Evo-POCS is classic POCS constrained by local dip consistency: it moves a missing sample toward a local wavefield prediction only where one dip direction is clearly more self-consistent than the others. Algorithm 2 summarizes the deployed Evo-POCS at the frequency-slice level. In each iteration, after the classic POCS sparse projection, a local phase relation is estimated from observed neighboring traces for a fixed set of candidate dip directions; each direction predicts the missing sample by phase-aligned continuation; the directions whose predictions best match the observed samples are combined into a consensus prediction; and a confidence weight, set by how clearly the best direction outperforms the rest, controls how far the estimate moves toward this prediction. The observed traces are restored exactly at each step. Implementation details of the phase coefficient, residual weighting, and fallback rules are given in Appendix A. Algorithm 2 Evo-POCS reconstruction for one frequency slice 1: observed slice YkY_k, mask M, iterations J, direction bank D 2: reconstructed slice UkU_k 3: estimate local phase coefficients cd(u)c_d(u) for each direction d∈d from observed trace pairs ⊳ depends only on YkY_k and M 4: Uk(0)←YkU_k^(0)← Y_k 5: for j=0,…,J−1j=0,…,J-1 do 6: V←V← Fourier hard-threshold sparse projection of Uk(j)U_k^(j) ⊳ classic POCS step 7: for each direction d∈d do 8: predict PdP_d from V by phase-aligned continuation 9: measure residual RdR_d on observed samples 10: end for 11: form consensus P⋆P_ from selected lowest-residual directions 12: set confidence map κ from residual contrast 13: V←V+κ⊙(P⋆−V)V← V+κ (P_ -V) ⊳ dip-consistency correction 14: Uk(j+1)←Yk+(1−M)⊙VU_k^(j+1)← Y_k+(1-M) V ⊳ restore observed traces 15: end for 16: return Uk(J)U_k^(J) After all frequency slices are reconstructed, an inverse temporal Fourier transform gives the real-valued volume. The correction is not an explicit isotropic smoothing operator, nor does it retune the threshold; instead, it uses the directional continuity present in the data itself, deciding along which direction and by how much to correct each missing sample directly from the observed data. Every part of the operator has a clear seismic meaning: the Fourier threshold corresponds to the sparsity assumption, the observed-trace projection keeps the measured data unchanged, the candidate direction set represents the local dip of seismic events, the phase estimate compensates the phase difference between neighboring traces, and the residual and confidence ensure that the dip prior is trusted only along the directions the observed data support. Unlike existing methods that rely on dip continuity, the operator in Evo-POCS was discovered autonomously by SeisEvo: without any prescribed mechanism, the search constructed the full procedure of directional estimation, consensus, and gating directly from the observed data, and integrated it into the POCS framework with a clear improvement. Reconstruction results. Figure 4 compares classic POCS and Evo-POCS on a 64×64×6464× 64× 64 sub-block under 50% random sampling. Panel (a) is the original complete data, (d) the 50% randomly sampled data, and (b) and (c) the reconstructions by classic POCS and Evo-POCS, respectively. Both methods reconstruct the data effectively, but Evo-POCS clearly yields better continuity of the local events, and is notably smoother on the time slices. The error panels (e) and (f) further highlight this difference: Evo-POCS produces less leakage, a smaller reconstruction error, and lower residual noise. Quantitatively, the SNR on this sub-block improves from 11.78 dB for classic POCS to 18.26 dB, SSIM from 0.895 to 0.986, and MSE decreases from 2.3×10−32.3× 10^-3 to 5.0×10−45.0× 10^-4. Figure 4: Reconstruction comparison on a 64×64×6464× 64× 64 synthetic sub-block under 50% random sampling. (a) Original complete data; (b) POCS reconstruction; (c) Evo-POCS reconstruction; (d) 50% randomly sampled data; (e) POCS error; (f) Evo-POCS error. Generalization test. To further assess the generalization of Evo-POCS, we test reconstruction on several out-of-distribution synthetic and field datasets. On these data, which were not used during the search, Evo-POCS improves over classic POCS by more than 3 dB on average. Figure 5 shows the reconstruction on one field dataset (256×120×64256× 120× 64) under 50% random sampling. Evo-POCS again gives clearly better event continuity and a smaller reconstruction error. The improvement is especially visible above the first arrivals, where Evo-POCS suppresses the noise introduced by classic POCS and markedly reduces the error in that region. Quantitatively, the SNR improves from 9.23 dB for classic POCS to 12.79 dB, and SSIM from 0.980 to 0.992. Figure 6 further shows the reconstruction on a slice at Y=0.65Y=0.65 km together with the corresponding F-K spectra. Both methods reconstruct the data effectively, but the advantage is clearer in this slice: Evo-POCS gives stronger event continuity and less noise, and its F-K spectrum is closer to that of the original data, with the sampling aliasing better suppressed. Figure 5: Reconstruction comparison on a 256×120×64256× 120× 64 field dataset under 50% random sampling. (a) Original; (b) POCS; (c) Evo-POCS; (d) sampled data; (e) POCS error; (f) Evo-POCS error. The Evo-POCS reconstruction error is smaller than that of classic POCS. Figure 6: Reconstruction and F-K comparison on a slice at Y=0.65Y=0.65 km. (a)–(d) original, sampled, POCS, and Evo-POCS reconstructions; (e)–(h) their F-K spectra. Evo-POCS gives better event continuity and an F-K spectrum closer to the original. We further evaluate the sensitivity to the missing ratio on the same field dataset. Table 1 reports the SNR comparison between POCS and Evo-POCS for missing ratios from 30% to 70%. Evo-POCS consistently outperforms POCS at all missing ratios, with SNR gains from 1.58 to 5.42 dB and an average gain of 3.49 dB. This result shows that the improvement is not limited to a single sampling density. The SNR gain gradually decreases as the missing ratio becomes more severe. This behavior is expected because the evolved dip-consistency correction relies on the structural evidence carried by the observed traces; when fewer traces are available, the local dip information becomes less reliable and the prior has weaker support. Nevertheless, Evo-POCS still improves over POCS even when 70% of the traces are missing. Table 1: SNR comparison between classic POCS and Evo-POCS on the field dataset under different missing-trace ratios. Δ is computed as SNREvo-POCS−SNRPOCSSNR_Evo -POCS-SNR_POCS. Missing ratio (%) POCS (dB) Evo-POCS (dB) Δ (dB) 30 9.62 15.04 + 5.42 40 9.48 13.73 + 4.25 50 9.23 12.79 + 3.56 60 8.88 11.52 + 2.64 70 8.41 9.99 + 1.58 Average 9.12 12.61 + 3.49 3.3 SeisEvo test on MSSA for simultaneous interpolation and denoising Task instantiation. The second case study instantiates SeisEvo on the MSSA family for simultaneous interpolation and denoising of 3D seismic data. The observations follow the noisy sampling model in Eq. (1), Y=M⊙(X+N)Y=M (X+N), where M is a random spatial sampling mask that retains 50% of the complete traces. The entries of N are independent zero-mean Gaussian variables, N∼(0,σ2)N (0,σ^2), with σ=αmax|X|σ=α |X| and α∼(0.01,0.20)α (0.01,0.20). Thus, the noise level is scaled to the amplitude of each data block and varies across the search cases. As in the POCS case study, the search uses 100 seismic sub-blocks of size 64×64×6464× 64× 64. Candidate algorithms are generated and scored only on these sub-blocks, and the external test data are not used during the search. The seed A0A_0 is a fixed-window implementation of the classical block-Hankel MSSA low-rank model. The editable surface S includes singular-value processing and rank selection, the Hankel low-rank projection, and the rule used to estimate observation reliability in the soft data update. The temporal Fourier transform, Hankel embedding, overlap averaging, input–output interface, and mask-based form of the data update remain fixed. The constraint set C requires the sampling mask to remain unchanged and the output to be finite and real-valued. It also forbids test-set leakage and unexplained free parameters. Unlike the POCS case without added synthetic noise, this joint task assigns different roles to the two data subsets: the acquired traces require denoising, whereas the missing traces require interpolation. We therefore separate the SNR-related assessment into the observed-trace denoising gain GobsG_obs and the missing-trace reconstruction SNR SNRmissSNR_miss, while SSIM measures structural fidelity over the complete reconstruction: OMSSA=0.35clip(Gobs18,0,1)+0.35clip(SNRmiss18,0,1)+0.30SSIM,O_MSSA=0.35\,clip ( G_obs18,0,1 )+0.35\,clip ( SNR_miss18,0,1 )+0.30\,SSIM, (12) where GobsG_obs is the SNR gain of the reconstructed observed traces over the noisy input, and SNRmissSNR_miss is evaluated only on the missing traces. The 18 dB reference scales the two SNR terms to ranges comparable with SSIM during the search; it is not a prescribed reconstruction target. MSSA seed. In practical seismic processing, MSSA is often applied within local spatial windows. Local processing limits the size of the block-Hankel matrix and reduces the cost of repeated low-rank decomposition. It also makes the low-rank approximation more appropriate because seismic events are closer to locally linear within a limited spatial aperture (31). Following this practice, we use a computationally compact, fixed-window implementation of classical block-Hankel MSSA (23) as the seed. Here, “compact” refers only to the use of a small local Hankel embedding. Throughout this paper, we denote this seed simply as MSSA. After a temporal Fourier transform, each frequency slice is processed independently. For the kkth frequency slice Uk(j)U_k^(j), its low-rank estimate is Lk(j)=ℋ†[(ℋ(Uk(j)))],L_k^(j)=H [P (H\! (U_k^(j) ) ) ], (13) where ℋH is the spatial block-Hankel embedding, P is the low-rank projection, and ℋ†H maps the projected Hankel matrix back to the spatial frequency slice by overlap averaging. The MSSA seed and all candidate algorithms use a fixed 6×66× 6 Hankel window, a rank limit of 16, and 30 reconstruction iterations. Because the observed traces are noisy, the seed does not restore them exactly. Instead, it applies the soft data-consistency update Uk(j+1)=M⊙[γk,jYk+(1−γk,j)Lk(j)]+(1−M)⊙Lk(j),U_k^(j+1)=M [ _k,jY_k+ (1- _k,j )L_k^(j) ]+(1-M) L_k^(j), (14) where γk,j∈[0,1] _k,j∈[0,1] balances the noisy observation YkY_k and the low-rank estimate Lk(j)L_k^(j). This implementation retains the central low-rank model of classical MSSA while allowing the acquired traces to be corrected during denoising. Search process. Starting from the MSSA seed, SeisEvo runs for 20 search rounds. Figure 7 shows the search-time objective OMSSAO_MSSA over the search. The score increases from 0.320 for the seed to 0.563 for the best candidate. Correspondingly, the observed-trace denoising gain GobsG_obs increases from 3.06 to 7.71 dB, the missing-trace reconstruction SNR SNRmissSNR_miss from 3.83 to 8.59 dB, and the overall SSIM from 0.619 to 0.819. Although the overall reconstruction SNR is not included directly in the search objective, it also increases from 4.37 to 9.17 dB. These concurrent gains show that the search improves observed-trace denoising, missing-trace interpolation, and structural fidelity together, rather than improving one component at the expense of another. Figure 7: Search-time objective OMSSAO_MSSA over 20 search rounds. The score increases from 0.320 for MSSA to 0.563 for the best candidate. Evo-MSSA. The main structural change found by SeisEvo is a reliability-grouped singular-value shrinkage. In one sentence, Evo-MSSA is MSSA with data-driven patch grouping: rather than fitting all Hankel patches with one fixed-rank subspace, it groups them by local observation reliability and fits each group with its own low-rank subspace. The fixed Hankel window, rank limit, and reconstruction count remain unchanged from the MSSA seed. Algorithm 3 summarizes the deployed operator for one frequency slice. The noisy observation YkY_k is fixed throughout the iterations and serves as both the initial value and the anchor of the soft data update. In the first iteration, when no previous reliability field is available, Evo-MSSA applies a global MP–GD shrinkage. Here, MP–GD denotes a singular-value shrinkage whose noise scale is estimated from the lower spectral tail and calibrated by the Marchenko–Pastur law, while its threshold coefficient is given by the Gavish–Donoho optimal hard-threshold form (19; 12). From the second iteration onward, the Hankel columns are divided at the empirical mean of the patch reliability derived from Γk(j−1) _k^(j-1). Each group is projected using its own singular spectrum and matrix aspect ratio, and the projected columns are placed back in their original positions. The resulting low-rank slice is then compared with YkY_k on the observed traces. Overlapping-support averaging and a centered cross-iteration Cesàro average convert these residuals into the current reliability field Γk(j) _k^(j), which controls the soft update. The complete MP–GD, grouping, reliability, and fallback expressions are given in Appendix B. Algorithm 3 Evo-MSSA reconstruction for one frequency slice 1: fixed noisy observed slice YkY_k; mask M; iterations J; Hankel operator ℋH 2: reconstructed slice UkU_k 3: Uk(0)←YkU_k^(0)← Y_k; Γk(−1)←∅ _k^(-1)← 4: for j=0,…,J−1j=0,…,J-1 do 5: πj←j/(J−1) _j← j/(J-1) 6: Bk(j)←ℋ(Uk(j))B_k^(j) (U_k^(j)) 7: if j=0j=0 then 8: B~k(j)←MPGD(Bk(j)) B_k^(j) (B_k^(j)) 9: else 10: B~k(j)←GroupMPGD(Bk(j),Γk(j−1)) B_k^(j) (B_k^(j), _k^(j-1)) 11: end if 12: Lk(j)←ℋ†(B~k(j))L_k^(j) ( B_k^(j)) 13: Γk(j)←ConsensusTrust(Yk,Lk(j),M,πj) _k^(j) (Y_k,L_k^(j),M, _j) 14: Uk(j+1)←M⊙[Γk(j)⊙Yk+(1−Γk(j))⊙Lk(j)]+(1−M)⊙Lk(j)U_k^(j+1)← M [ _k^(j) Y_k+(1- _k^(j)) L_k^(j)]+(1-M) L_k^(j) 15: end for 16: return Uk(J)U_k^(J) The three modules in Algorithm 3 have distinct roles. MPGDMPGD adapts the shrinkage to the measured spectrum and matrix shape; GroupMPGDGroupMPGD permits patches with different reliability to use different subspaces; and ConsensusTrustConsensusTrust stabilizes the residual-based reliability estimate in space and across iterations. The evolved mechanism introduces no manually tuned continuous parameter. The two-group structure, empirical-mean split, robust residual scale, and MP plug-in remain explicit structural choices rather than assumptions hidden in learned weights. Reconstruction results. Figure 8 compares MSSA and Evo-MSSA on a 64×64×6464× 64× 64 synthetic sub-block with 50% randomly missing traces and Gaussian noise. Panel (a) shows the original complete data, panel (d) the noisy incomplete data, panels (b) and (c) the reconstructions obtained by MSSA and Evo-MSSA, and panels (e) and (f) their corresponding errors. Both methods recover the main seismic events across the missing traces. However, visible residual noise remains in the MSSA result. Evo-MSSA produces more continuous events, weaker background noise, and a smoother reconstruction. The error panels further show that Evo-MSSA leaves less residual noise and a smaller overall reconstruction error. The input SNR of the noisy incomplete data is 1.91 dB. MSSA increases the SNR to 7.03 dB, whereas Evo-MSSA reaches 11.09 dB, an improvement of 4.06 dB over MSSA. The MSE is also reduced from 6.74×10−36.74× 10^-3 for MSSA to 2.65×10−32.65× 10^-3 for Evo-MSSA. These results show that the Evo-MSSA operator obtained through the SeisEvo search clearly improves on the MSSA seed. It recovers events across missing traces while more effectively attenuating noise on the observed traces. Figure 8: Reconstruction comparison on a 64×64×6464× 64× 64 synthetic sub-block with 50% randomly missing traces and Gaussian noise. (a) Original complete data; (b) MSSA reconstruction; (c) Evo-MSSA reconstruction; (d) noisy data with 50% randomly missing traces; (e) MSSA error; (f) Evo-MSSA error. Generalization test. To further assess the generalization of Evo-MSSA, we test simultaneous interpolation and denoising on several synthetic and field datasets that were not used during the search. Across these unseen datasets, Evo-MSSA improves the overall reconstruction SNR over MSSA by more than 4 dB on average. This result shows that the mechanism found by SeisEvo is not limited to the sub-blocks used during the search and can generalize to unseen synthetic and field seismic data. Figure 9 shows a representative result on a 256×120×64256× 120× 64 field dataset outside the search set. Panel (a) shows the complete data after noise injection, whereas panel (d) shows the same noisy data after 50% of the traces are removed. The SNR of the noisy incomplete input is 2.95 dB. MSSA increases the reconstruction SNR to 10.88 dB, whereas Evo-MSSA reaches 14.49 dB, a gain of 3.61 dB over MSSA. The corresponding MSE decreases from 9.84×10−59.84× 10^-5 for MSSA to 3.92×10−53.92× 10^-5 for Evo-MSSA. Both methods recover the main seismic events, but the Evo-MSSA result is smoother and shows stronger noise attenuation. The difference is clearer in the error panels: visible residual noise remains after MSSA reconstruction, whereas Evo-MSSA produces a weaker and cleaner reconstruction error. Figure 9: Reconstruction comparison on a 256×120×64256× 120× 64 field dataset with 50% randomly missing traces and Gaussian noise. (a) Noisy complete data; (b) MSSA reconstruction; (c) Evo-MSSA reconstruction; (d) noisy incomplete data; (e) MSSA error; (f) Evo-MSSA error. Comparison with ODRR. To determine whether the improvement extends beyond the seed method, we further compare Evo-MSSA with the optimally damped rank-reduction method (ODRR) (5; 6). Following the recommended parameter range, a small validation sweep selected K=3K=3. Except for the damping setting specific to ODRR, the reconstruction parameters shared by the methods were kept unchanged. We vary the noise scale α from 0.01 to 0.05 on the same field dataset with 50% randomly missing traces. Table 2 reports the overall reconstruction SNR. Evo-MSSA achieves the highest SNR at every noise level, and its advantage becomes more pronounced at higher noise levels. Averaged over the five noise levels, Evo-MSSA reaches 10.70 dB, exceeding MSSA and ODRR by 7.15 and 3.31 dB, respectively. These results show that SeisEvo improves the classical MSSA seed and remains effective relative to the ODRR baseline. Table 2: Overall reconstruction SNR on the field dataset under different noise scales. Noise scale α MSSA (dB) ODRR (dB) Evo-MSSA (dB) Gain over ODRR (dB) 0.01 10.88 12.61 14.49 +1.88+1.88 0.02 5.77 8.84 12.05 +3.21+3.21 0.03 2.58 7.13 10.32 +3.19+3.19 0.04 0.23 5.50 8.96 +3.46+3.46 0.05 -1.71 2.85 7.70 +4.85+4.85 Average 3.55 7.39 10.70 +3.31+3.31 Figure 10 provides a visual comparison at α=0.05α=0.05. Under this low-input-SNR condition, simultaneous interpolation and denoising are difficult. The MSSA reconstruction remains dominated by noise, and most seismic events are barely identifiable. ODRR recovers part of the event structure, but substantial residual noise and spatial artifacts remain. Evo-MSSA gives the clearest result: the main events remain visible and show better continuity. Nevertheless, its error volume still contains coherent signal energy and random noise, indicating signal leakage and residual noise. This comparison shows that Evo-MSSA is more robust than MSSA and ODRR under this low-input-SNR condition, while also revealing its limitations in this challenging case. Figure 10: Reconstruction comparison on a 256×120×64256× 120× 64 field dataset at noise scale α=0.05α=0.05 with 50% randomly missing traces. (a) Noisy complete data; (b) MSSA reconstruction; (c) ODRR reconstruction; (d) Evo-MSSA reconstruction; (e) noisy sampled data; (f) MSSA error; (g) ODRR error; (h) Evo-MSSA error. Figure 11 further compares the slice at Y=1.75Y=1.75 km and the corresponding F–K spectra. At α=0.05α=0.05, random noise and missing traces make the original events almost unidentifiable in the noisy sampled data. Its F–K spectrum also contains little identifiable signal energy. MSSA does not effectively separate the signal from noise: its reconstruction remains dominated by noise, and the corresponding F–K spectrum shows the same behavior. ODRR recovers part of the main events and their spectral structures, but substantial residual noise remains. In comparison, Evo-MSSA recovers clearer and more continuous main events and further suppresses the effects of random noise and missing traces. Its F–K energy is more concentrated around the main signal structures and is closer to that of the original data. However, some weak events and spectral energy are still attenuated, indicating remaining signal loss under this low-input-SNR condition. The slice comparison is consistent with the full-volume results and quantitative metrics, further showing the improved interpolation and denoising performance of Evo-MSSA at low SNR. Figure 11: Slice and F–K spectrum comparison on a 256×120×64256× 120× 64 field dataset at Y=1.75Y=1.75 km, with noise scale α=0.05α=0.05 and 50% randomly missing traces. (a) Original; (b) noisy sampled data; (c) MSSA; (d) ODRR; (e) Evo-MSSA; (f)–(j) the corresponding F–K spectra of (a)–(e), respectively. Taken together, the two case studies start from Fourier-sparse POCS and low-rank MSSA and yield Evo-POCS and Evo-MSSA, respectively. In the POCS case, Evo-POCS retains the Fourier sparse projection and data-consistency constraint while mainly adding a data-driven dip-consistency correction, which improves the continuity of reconstructed events. In the MSSA case, Evo-MSSA improves the conventional low-rank reconstruction through adaptive singular-value shrinkage, reliability grouping, and cross-iteration consensus. It also maintains an advantage over ODRR, especially at higher noise levels. The two cases address different reconstruction tasks: Evo-POCS considers interpolation of randomly missing traces without added synthetic noise, whereas Evo-MSSA considers simultaneous interpolation and denoising with randomly missing traces and Gaussian noise. In their respective tasks, both evolved algorithms clearly outperform their classical seeds and maintain consistent improvements on data outside the search set. These results indicate that the mechanisms discovered by SeisEvo can work with algorithm families based on different mathematical principles and under different observation conditions and reconstruction objectives, rather than being limited to the sub-blocks used during search. Overall, the experiments support the effectiveness of SeisEvo on two classical seismic reconstruction families and two different reconstruction tasks. 4 Discussion What the search actually contributed. The discovered mechanisms are not conceptually new. Interpolation that exploits dip continuity, and adaptive shrinkage of singular values, each have their own research background in seismic processing. The value of SeisEvo does not lie in proposing a previously unseen prior, but in the fact that the specific form of these mechanisms was constructed by the search without the form of the operator being prescribed, and was integrated automatically into the iteration of a classical solver. In Evo-POCS, the directional candidate set, the phase estimation from observed trace pairs, the residual-weighted directional consensus, and the confidence gating were not given by the specification; in Evo-MSSA, neither was the structure that divides the Hankel columns by local reliability and projects each group with its own subspace. These elements were built up and selected under the physical-legality and parameter-provenance constraints. The resulting operators remain explicit algorithms that can be written down, inspected at the operator level, related to interpretable seismic mechanisms, and tested independently. This is what distinguishes algorithm evolution from replacing a solver with an implicit mapping: it changes the way an algorithm is designed, not the way an algorithm is expressed. Scope and cost of the protocol. More generally, the conditions on which SeisEvo relies, namely an executable quality measure, explicit physical contracts, and a classical algorithm with a clear and modular structure, are also present in seismic denoising, deconvolution, and inversion, so the protocol may extend to other steps of seismic signal processing. The main cost at present is offline search computation: discovering an operator requires many candidate executions and is therefore substantially more expensive than running a fixed hand-designed solver. This cost is paid only during discovery; once frozen, the operator requires no agent, prompt, or neural network at inference time. The present evidence, however, is restricted to seismic data reconstruction; extension to these tasks remains a hypothesis to be tested rather than a demonstrated capability. Design of the objective and evaluator. This work makes the constraint set and the editable operator surface explicit and part of the search protocol, but the seed and the objective are fixed before each run, and the effect of their choice has not been studied systematically. The objective encodes part of the human judgement of quality as an executable scoring rule. The two case studies use different scoring structures: the POCS search combines the reconstruction SNR with the structural similarity, whereas the MSSA search further separates noise suppression at the observed positions from reconstruction quality at the missing ones. The associated weights and normalizing references belong to an evaluation protocol that is identical for all candidates; they are fixed before the search and held by the evaluator, and no candidate program may modify them. The objective plays a role comparable to the loss function of a learning method, in that both prescribe the optimization preference, but it guides a search over discrete programs and operator mechanisms rather than the optimization of network weights. It therefore need not be differentiable and may include piecewise measures and discrete quality judgements, while legality checks are enforced separately through the constraint set. Different objectives consequently create different search pressures: raising the SNR, preserving structure, suppressing noise, and limiting signal leakage are not always mutually consistent. An executable evaluator thus does more than record performance; it defines the selection pressure under which a program evolves (22). The executable form of the objective can itself become an object of search, as in the automatic design and evolution of reward code (18). In the present runs, the best-so-far scores exhibit late-stage plateaus, but a single trajectory per case is not sufficient to decide whether the budget is already adequate, whether candidate diversity is being exhausted, or whether the room for observable improvement under the current evaluation protocol has simply become small. Separating these possibilities would require repeated independent searches, controlled experiments with longer budgets, and explicit stopping rules based on performance gain, candidate diversity, and computational cost. From offline discovery to continual discovery. The SNR and the structural similarity used here both require a complete reference, so the present protocol suits offline algorithm discovery when a complete reference is available or when degraded observations can be constructed from complete data, whereas the traces that are genuinely missing in an acquired survey are themselves unobservable. For interpolation, part of the observed traces can be hidden and used for self-supervised validation; for simultaneous interpolation and denoising, the hidden traces are themselves noisy, so additional assumptions about noise independence, repeated observations, or noise statistics are needed. Spectral consistency, local coherence, and the quality of subsequent imaging can also serve as substitutes, but they remain proxy objectives and may introduce their own bias into the search, which is consistent with the preceding observation that the scoring design shapes the mechanism that is found. A separate issue is the dependence on the seed and the search backend. Each case study here examines a single classical seed with a fixed backend; whether different seeds would arrive at comparable mechanisms, and how far objectives and constraints with more specific geophysical meaning, such as band preservation, amplitude fidelity, or a forward-modelling residual, would shape the direction of the search, remain to be examined. If the discovery process is to continue, for instance by taking a discovered operator as the seed of a further run, then repeatedly evaluating against the same external data would gradually make that data part of the search and cost it its role as independent validation. Beyond this, what evolves in this work is the reconstruction solver rather than the search mechanism itself. Adaptive optimization of candidate generation, experience organization, and budget allocation is a natural extension: as long as the evaluator and the hard constraints stay fixed, the legality decisions and the scores of executed candidates remain reproducible under fixed programmatic rules. If the objective or the legality rules were to co-evolve as well, however, the standard of evaluation would become a moving target, and a fixed outer validation protocol that does not participate in the evolution would be required; otherwise an increase in score would no longer necessarily correspond to an improvement of the algorithm itself. 5 Conclusions We proposed SeisEvo (Seismic Algorithm Evolution), an algorithm evolution framework for seismic data reconstruction. Rather than optimizing a single reconstructed volume, SeisEvo searches for the reconstruction algorithm itself. A domain expert defines a classical seed algorithm, the components that may be edited, the reconstruction objective, and the legality constraints; an LLM-based multi-agent search then improves executable candidates within those bounds. The final output is a standalone white-box operator that requires no agent, prompt, or neural network at inference time. We instantiated SeisEvo on two classical algorithm families with different mathematical foundations. For interpolation of field data without added synthetic noise, Evo-POCS improves the SNR over classic POCS by 3.49 dB on average across missing ratios from 30% to 70%, with the gain decreasing from 5.42 dB at 30% missing to 1.58 dB at 70% missing. For simultaneous interpolation and denoising, Evo-MSSA improves the average reconstruction SNR by more than 7 dB over classic MSSA and by more than 3 dB over the ODRR baseline, under 50% randomly missing traces and across the tested noise levels. Both evolved operators maintain their performance gains on synthetic and field data not used during the search. These results suggest a different role for domain expertise in geophysical algorithm design: the expert defines what constitutes a legal and useful algorithm, while the mechanism itself is discovered through search, and the resulting operators remain explicit, inspectable, and directly deployable. SeisEvo therefore indicates the potential of agentic algorithm evolution as a complementary path to deep learning for discovering interpretable geophysical processing algorithms. Acknowledgments This work was supported in part by Deep Earth Probe and Mineral Resources Exploration - National Science and Technology Major Project under Grant No. 2024ZD1002700, in part by the NSFC under Grant No. 42574172. Yingjie Xu acknowledges the support of the China Scholarship Council program (CSC No. 202506120138). Appendix A Mathematical formulation of the dip-consistency correction in Evo-POCS This appendix gives the mathematical form of the dip-consistency correction in Evo-POCS. All operations are performed on a single frequency slice. Let YkY_k be the observed frequency slice used in Eq. (11), M the sampling mask, U the current iterate, and u=(x,y)u=(x,y) a spatial location. Shifted locations u±du± d that fall outside the slice are excluded from the local sums below. Direction bank The correction is applied over a fixed set of candidate dip directions, each a small integer spatial offset: =(1,0),(0,1),(1,1),(1,−1),(2,1),(1,2),(2,−1),(1,−2),(2,0),(0,2),(2,2),(2,−2). splitD=\&(1,0),(0,1),(1,1),(1,-1),(2,1),(1,2),\\ &(2,-1),(1,-2),(2,0),(0,2),(2,2),(2,-2)\. split (A.1) These twelve offsets, with steps at most two, form a discrete set of local dip directions for seismic events. Local phase coefficient For each direction d=(dx,dy)d=(d_x,d_y), a complex phase coefficient cd(u)c_d(u) is estimated from observed trace pairs, describing the phase rotation between neighboring traces along that direction. Within a 7×77× 7 window 7(u)N_7(u) centered at u, cd(u)=∑v∈7(u)Yk(v+d)¯Yk(v)M(v+d)M(v)|∑v∈7(u)Yk(v+d)¯Yk(v)M(v+d)M(v)|+ϵ,c_d(u)= _v _7(u) Y_k(v+d)\,Y_k(v)\,M(v+d)M(v) | _v _7(u) Y_k(v+d)\,Y_k(v)\,M(v+d)M(v) |+ε, (A.2) where ϵε is a small constant that prevents division by zero. The normalization makes cd(u)c_d(u) have approximately unit modulus, so it acts mainly as a phase rotation; the mask product M(v+d)M(v)M(v+d)M(v) ensures that only trace pairs observed at both ends are used. Since cd(u)c_d(u) depends only on the observed slice and the mask, it is computed once before the iterations. Directional prediction and residual Given the current iterate U, each direction predicts the current location by symmetric phase-aligned continuation from its two neighbors: Pd()(u)=12[cd(u)(u+d)+cd(u)¯(u−d)].P_d(U)(u)= 12 [c_d(u)\,U(u+d)+ c_d(u)\,U(u-d) ]. (A.3) Indeed, if the local wavefield satisfies (u+d)=eiϕd(u)U(u+d)=e^i _dU(u), then cd(u)≈e−iϕdc_d(u)≈ e^-i _d under Eq. (A.2). Consequently, both cd(u)(u+d)c_d(u)U(u+d) and cd(u)¯(u−d) c_d(u)U(u-d) predict (u)U(u); the conjugate in the second term provides the reverse-direction phase alignment. The reliability of this prediction is measured by its residual on observed samples, within a 3×33× 3 window 3(u)N_3(u): Rd(u)=∑v∈3(u)M(v)|Yk(v)−Pd()(v)|2∑v∈3(u)M(v)+ϵ.R_d(u)= _v _3(u)M(v)\,|Y_k(v)-P_d(U)(v)|^2 _v _3(u)M(v)+ε. (A.4) A smaller residual indicates that the direction better matches the local dip of the events. Consensus prediction Let R¯(u) R(u) be the mean residual over the twelve directions, and let d1,d2,d3d_1,d_2,d_3 be the three directions with the smallest residuals. They are combined into a consensus prediction using cubic residual-margin weights: wi(u)=[R¯(u)−Rdi(u)]+3,P⋆()(u)=∑i=13wi(u)Pdi()(u)∑i=13wi(u)+ϵ,w_i(u)= [ R(u)-R_d_i(u) ]_+^3, P_ (U)(u)= _i=1^3w_i(u)\,P_d_i(U)(u) _i=1^3w_i(u)+ε, (A.5) where [⋅]+[·]_+ takes the non-negative part. Directions with smaller residuals receive larger weights, so the most self-consistent direction dominates the prediction. Confidence gate and update The strength of the correction is controlled by a confidence κ(u)κ(u), which measures how much the best direction improves over the mean residual: κ(u)=clip(R¯(u)−Rmin(u)R¯(u)+ϵ, 0, 1),Rmin(u)=mind∈Rd(u).κ(u)=clip\! ( R(u)-R_ (u) R(u)+ε,\,0,\,1 ), R_ (u)= _d R_d(u). (A.6) When one direction clearly outperforms the rest, κ is close to one and the correction is strong; when the residuals are similar, κ is close to zero and the correction is weak. Writing POCSU^POCS for the output of the classic POCS sparse projection (i.e., Vk(j)V_k^(j) in Eq. (10)), one full Evo-POCS update is (j+1)=Yk+(1−M)⊙[POCS+κ⊙(P⋆(POCS)−POCS)].U^(j+1)=Y_k+(1-M) [U^POCS+κ (P_ (U^POCS)-U^POCS ) ]. (A.7) The outer term Yk+(1−M)⊙(⋅)Y_k+(1-M) (·) restores the observed traces exactly, so only the missing traces are updated by the correction. Fallback rules When a 7×77× 7 window contains no usable observed trace pair, cd(u)c_d(u) falls back to a global phase estimate for the slice; when a 3×33× 3 window contains no observed sample, the residual Rd(u)R_d(u) is replaced by a self-consistency measure across directions. These fallbacks help keep the operator stable where observations are sparse. Real output After all slices are processed, Hermitian symmetry is imposed before the inverse temporal Fourier transform, which gives a real-valued volume. Appendix B Mathematical formulation of Evo-MSSA This appendix gives the detailed form of the three modules summarized in Algorithm 3. All operations below are performed on one temporal frequency slice. The noisy observed slice YkY_k and mask M remain fixed over the 30 reconstruction iterations, whereas Uk(j)U_k^(j) is the current iterate. Local Hankel embedding For a slice Uk(j)U_k^(j), the operator ℋH extracts all overlapping 6×66× 6 spatial patches, vectorizes them, and places them as columns of Bk(j)=ℋ(Uk(j))=[b1,…,bnp],bp∈ℂ36.B_k^(j)=H(U_k^(j))= [b_1,…,b_n_p ], b_p ^36. (B.1) Here, npn_p is the number of overlapping patches. The associated operator ℋ†H returns the columns to their original patch locations and forms a weighted overlap-add reconstruction using a separable Hanning taper on the 6×66× 6 support, normalized pointwise by the accumulated taper weights. Thus, ℋ†H denotes tapered overlap-averaging reconstruction, not a Moore–Penrose pseudoinverse. MP–GD singular-value shrinkage Consider either the complete Hankel matrix or one of its column groups, denoted by Bg∈ℂmg×ngB_g ^m_g× n_g. Its left Gram decomposition is BgBgH=Qgdiag(sg,12,…,sg,rg2)QgH,sg,1≥⋯≥sg,rg≥0,B_gB_g^H=Q_gdiag (s_g,1^2,…,s_g,r_g^2 )Q_g^H, s_g,1≥·s≥ s_g,r_g≥ 0, (B.2) where (⋅)H(·)^H is the conjugate transpose and rg=min(mg,ng)r_g= (m_g,n_g). Let gT_g be the lower half of this singular spectrum, and let ηg=min(mg,ng)/max(mg,ng) _g= (m_g,n_g)/ (m_g,n_g) be the measured matrix aspect ratio. If qgMPq_g^MP denotes the 25th percentile of the Marchenko–Pastur eigenvalue law (19), the matrix noise scale is σ^B,g=median(g)qgMP. σ_B,g= median(T_g) q_g^MP. (B.3) The 25th percentile appears because the median of the lower half of the full spectrum is its first quartile. If the plug-in estimate is unavailable, the implementation uses the largest finite value among the tail median, its MAD-based scale, and the numerical floor ϵε. Let ω⋆(ηg) _ ( _g) be the closed-form Gavish–Donoho optimal hard-threshold coefficient (12). Evo-MSSA sets θg=ω⋆(ηg)σ^B,g _g= _ ( _g) σ_B,g and forms the quadratic shrinkage weights wg,i=[1−(θgsg,i+ϵ)2]+.w_g,i= [1- ( _gs_g,i+ε )^2 ]_+. (B.4) Let rg⋆=min(16,rg)r_g = (16,r_g) be the inherited rank limit after accounting for the available matrix rank, and let Qg,rg⋆Q_g,r_g contain the corresponding leading left singular vectors. The MP–GD projection is MPGD(Bg)=Qg,rg⋆diag(wg,1,…,wg,rg⋆)Qg,rg⋆HBg,MPGD(B_g)=Q_g,r_g \,diag(w_g,1,…,w_g,r_g )\,Q_g,r_g ^HB_g, (B.5) with components beyond the inherited rank limit assigned zero weight. In the implementation, trailing zero-weight components are omitted. If all weights vanish, the projected group is zero, after which overlap averaging and the soft update proceed unchanged. The MP law therefore calibrates the noise scale, the Gavish–Donoho form sets the threshold location, and Eq. (B.4) performs the actual shrinkage. These quantities are used only as data-dependent calibration scales; no optimality is claimed for the resulting quadratic shrinkage on the correlated Hankel matrices considered here. Reliability-grouped projection No previous reliability field is available at j=0j=0, so Eq. (B.5) is applied to the complete matrix Bk(0)B_k^(0). For j≥1j≥ 1, let pW_p be the spatial support of column bpb_p. Its mean reliability and the resulting two groups are ξp(j−1) _p^(j-1) =1|p|∑u∈pΓk(j−1)(u), = 1|W_p| _u _p _k^(j-1)(u), ξ¯(j−1) ξ^(j-1) =1np∑p=1npξp(j−1), = 1n_p _p=1^n_p _p^(j-1), (B.6) ℐ−(j) _-^(j) =p:ξp(j−1)<ξ¯(j−1), =\p: _p^(j-1)< ξ^(j-1)\, ℐ+(j) _+^(j) =p:ξp(j−1)≥ξ¯(j−1). =\p: _p^(j-1)≥ ξ^(j-1)\. Each group is projected using its own spectrum and aspect ratio, and its columns are then returned to their original positions: B~k,:,ℐ±(j)(j)=MPGD(Bk,:,ℐ±(j)(j)). B_k,:,\,I_±^(j)^(j)=MPGD (B_k,:,\,I_±^(j)^(j) ). (B.7) This operation is denoted by GroupMPGDGroupMPGD in Algorithm 3. Residual-based reliability Let Ek(j)=Yk−Lk(j)E_k^(j)=Y_k-L_k^(j) be the prediction residual, evaluated only where M=1M=1. Let K be the 6×66× 6 all-ones support kernel and ∗* denote reflected-boundary spatial convolution. The local observed mass, residual power, and observed power are mK=K∗M,PE=K∗(M⊙|Ek(j)|2),PY=K∗(M⊙|Yk|2).m_K=K*M, P_E=K* (M |E_k^(j)|^2 ), P_Y=K* (M |Y_k|^2 ). (B.8) The global MP noise scale from the first projection is retained as a per-frequency anchor. With dM=(nxny)−1∑uM(u)d_M=(n_xn_y)^-1 _uM(u) denoting the sampling density, the corresponding noise power on observed samples is v^k=(σ^B,k(0))2npdM+ϵ. v_k= ( σ_B,k^(0))^2n_p\,d_M+ε. (B.9) Two bounded local signal fractions are then aobs=clip(1−PEPY+ϵ,0,1),anoise=clip(1−v^kPE/(mK+ϵ)+ϵ,0,1).a_obs=clip (1- P_EP_Y+ε,0,1 ), a_noise=clip (1- v_kP_E/(m_K+ε)+ε,0,1 ). (B.10) Writing πj=j/(J−1) _j=j/(J-1) for iteration progress, the raw local reliability and its overlapping-support average are araw(j)=(1−πj)aobs+πjanoise,asp(j)=K∗(M⊙araw(j))mK+ϵ.a_raw^(j)=(1- _j)a_obs+ _ja_noise, a_sp^(j)= K*(M a_raw^(j))m_K+ε. (B.11) Cross-iteration consensus Let the observed-sample mean and centered spatial deviation be μj=∑uM(u)asp(j)(u)∑uM(u)+ϵ,Dj=asp(j)−μj. _j= _uM(u)a_sp^(j)(u) _uM(u)+ε, D_j=a_sp^(j)- _j. (B.12) Evo-MSSA forms a centered Cesàro average across iterations, D¯j=D0,j=0,D¯j−1+Dj−D¯j−1j+1,j≥1,Γloc(j)=clip(μj+D¯j,0,1). D_j= casesD_0,&j=0,\\ D_j-1+ D_j- D_j-1j+1,&j≥ 1, cases _loc^(j)=clip( _j+ D_j,0,1). (B.13) This average reduces iteration-to-iteration fluctuations without introducing a forgetting factor. For robustness to complex residual outliers, let ℰk(j)=Ek(j)(u):M(u)=1E_k^(j)=\E_k^(j)(u):M(u)=1\, nE=|ℰk(j)|n_E=|E_k^(j)|, and cE=median(ℜℰk(j))+imedian(ℑℰk(j))c_E=median( _k^(j))+i\,median( _k^(j)). For nE>0n_E>0, define vmean=nE−1∑e∈ℰk(j)|e−cE|2v_mean=n_E^-1 _e _k^(j)|e-c_E|^2 and vrob v_rob =median(|ℰk(j)−cE|2)log2, = median (|E_k^(j)-c_E|^2 ) 2, π~j π_j =πj+1−πjnE, = _j+ 1- _jn_E, (B.14) Γrob(j) _rob^(j) =(1−π~j)+π~jclip[1−(1−1nE)vrobvmean+ϵ,0,1]. =(1- π_j)+ π_j\,clip [1- (1- 1n_E ) v_robv_mean+ε,0,1 ]. Both πj _j and its finite-sample correction π~j π_j gradually move the update from conservative retention of observed samples to residual- and noise-calibrated weighting. The MP-anchored signal fraction and the final map used by ConsensusTrustConsensusTrust are ζk _k =clip(1−v^k∑uM(u)∑uM(u)|Yk(u)|2+ϵ,0,1), =clip (1- v_k _uM(u) _uM(u)|Y_k(u)|^2+ε,0,1 ), (B.15) Γk(j) _k^(j) =(1−ζk)Γloc(j)+ζkΓrob(j). =(1- _k) _loc^(j)+ _k _rob^(j). Soft update and fallback rules Overlap averaging gives Lk(j)=ℋ†(B~k(j))L_k^(j)=H ( B_k^(j)), and the next iterate is Uk(j+1)=M⊙[Γk(j)⊙Yk+(1−Γk(j))⊙Lk(j)]+(1−M)⊙Lk(j).U_k^(j+1)=M [ _k^(j) Y_k+(1- _k^(j)) L_k^(j) ]+(1-M) L_k^(j). (B.16) Thus, large reliability values retain more of the noisy observation, whereas small values place more weight on the low-rank prediction. The observed traces are never hard-restored. The first iteration, a malformed reliability field, nonfinite patch scores, a group with fewer columns than Hankel rows, or a failed group eigendecomposition triggers the global MP–GD projection. Invalid spectral estimates use the finite tail-median/MAD fallback. If the observed-residual set is empty, the slice-level reliability falls back to the observed signal fraction; an invalid reconstructed frequency slice falls back to the observed slice. The final operator uses all 30 reconstruction iterations. After all positive frequency slices are processed, the inverse real Fourier transform produces a real-valued reconstructed volume. References Abedi and Pardo (2022) M. M. Abedi and D. Pardo A multidirectional deep neural network for self-supervised reconstruction of seismic data. IEEE Transactions on Geoscience and Remote Sensing 60, p. 1–9. External Links: Document Cited by: §1. Abma and Kabir (2006) R. Abma and N. Kabir 3D interpolation of irregular data with a POCS algorithm. Geophysics 71 (6), p. E91–E97. External Links: Document Cited by: §1, §1, §3.2. Bayati and Trad (2023) F. Bayati and D. Trad 3-D data interpolation and denoising by an adaptive weighting rank-reduction method using multichannel singular spectrum analysis algorithm. Sensors 23 (2), p. 577. Cited by: §2.1. Chen and Liu (2024) G. Chen and Y. Liu Combining unsupervised deep learning and Monte Carlo dropout for seismic data reconstruction and its uncertainty quantification. Geophysics 89 (1), p. WA53–WA65. External Links: Document Cited by: §1. Chen et al. (2020) Y. Chen, M. Bai, Z. Guan, Q. Zhang, M. Zhang, and H. Wang Five-dimensional seismic data reconstruction using the optimally damped rank-reduction method. Geophysical Journal International 222 (3), p. 1824–1845. External Links: Document Cited by: §1, §3.3. Chen et al. (2023) Y. Chen, W. Huang, L. Yang, Y. A. S. I. Oboué, O. M. Saad, and Y. Chen DRR: an open-source multi-platform package for the damped rank-reduction method and its applications in seismology. Computers & Geosciences 180, p. 105440. External Links: Document Cited by: §1, §3.3. Cheng et al. (2025a) M. Cheng, J. Lin, X. Dong, and T. Zhong A multitask deep-learning model for the denoising, interpolation, and wavefield separation of distributed acoustic sensing-vertical seismic profiling data. Geophysics 90 (6), p. V559–V568. External Links: Document Cited by: §1. Cheng et al. (2025b) S. Cheng, R. Harsuko, and T. Alkhalifah A generative foundation model for an all-in-one seismic processing framework. Surveys in Geophysics 46 (6), p. 1173–1215. External Links: Document Cited by: §1. Dong et al. (2025) L. Dong, M. Zhang, Y. Xu, C. Wang, S. Yu, and Y. Zhang Robust reconstruction of non-uniformly sampled 3D seismic data with outliers. Journal of Geophysics and Engineering 22 (2), p. 560–573. Cited by: §2.1. Gao et al. (2010) J. Gao, X. Chen, J. Li, G. Liu, and J. Ma Irregular seismic data reconstruction based on exponential threshold model of POCS method. Applied Geophysics 7 (3), p. 229–238. Cited by: §1, §1, §2.1, §3.2. Gao et al. (2013) J. Gao, A. Stanton, M. Naghizadeh, M. D. Sacchi, and X. Chen Convergence improvement and noise attenuation considerations for beyond alias projection onto convex sets reconstruction. Geophysical Prospecting 61 (S1), p. 138–151. External Links: Document Cited by: §1. Gavish and Donoho (2014) M. Gavish and D. L. Donoho The optimal hard threshold for singular values is 4/34/ 3. IEEE Transactions on Information Theory 60 (8), p. 5040–5053. External Links: Document Cited by: Appendix B, §3.3. Herrmann and Hennenfent (2008) F. J. Herrmann and G. Hennenfent Non-parametric seismic data recovery with curvelet frames. Geophysical Journal International 173 (1), p. 233–248. External Links: Document Cited by: §1. Huang et al. (2016) W. Huang, R. Wang, Y. Chen, H. Li, and S. Gan Damped multichannel singular spectrum analysis for 3D random noise attenuation. Geophysics 81 (4), p. V261–V270. Cited by: §1, §2.1. Kaur et al. (2021) H. Kaur, N. Pham, and S. Fomel Seismic data interpolation using deep learning with generative adversarial networks. Geophysical Prospecting 69 (2), p. 307–326. External Links: Document Cited by: §1. Liu et al. (2024) F. Liu, X. Tong, M. Yuan, X. Lin, F. Luo, Z. Wang, Z. Lu, and Q. Zhang Evolution of heuristics: towards efficient automatic algorithm design using large language model. In Proceedings of the 41st International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 235, p. 32201–32223. External Links: Link Cited by: §1. Liu et al. (2021) Q. Liu, L. Fu, and M. Zhang Deep-seismic-prior-based reconstruction of seismic data using convolutional neural networks. Geophysics 86 (2), p. V131–V142. External Links: Document Cited by: §1. Ma et al. (2024) Y. J. Ma, W. Liang, G. Wang, D. Huang, O. Bastani, D. Jayaraman, Y. Zhu, L. Fan, and A. Anandkumar Eureka: human-level reward design via coding large language models. In The Twelfth International Conference on Learning Representations, External Links: Link Cited by: §1, §4. Marchenko and Pastur (1967) V. A. Marchenko and L. A. Pastur Distribution of eigenvalues for some sets of random matrices. Mathematics of the USSR-Sbornik 1 (4), p. 457–483. External Links: Document Cited by: Appendix B, §3.3. Meng et al. (2022) F. Meng, Q. Fan, and Y. Li Self-supervised learning for seismic data reconstruction and denoising. IEEE Geoscience and Remote Sensing Letters 19, p. 1–5. External Links: Document Cited by: §1. Naghizadeh and Sacchi (2007) M. Naghizadeh and M. D. Sacchi Multistep autoregressive reconstruction of seismic records. Geophysics 72 (6), p. V111–V118. Cited by: §1, §1. Novikov et al. (2025) A. Novikov, N. Vũ, M. Eisenberger, E. Dupont, P. Huang, A. Z. Wagner, S. Shirobokov, B. Kozlovskii, F. J. Ruiz, A. Mehrabian, et al. AlphaEvolve: a coding agent for scientific and algorithmic discovery. arXiv preprint arXiv:2506.13131. External Links: Document Cited by: §1, §2.3, §4. Oropeza and Sacchi (2011) V. Oropeza and M. D. Sacchi Simultaneous seismic data denoising and reconstruction via multichannel singular spectrum analysis. Geophysics 76 (3), p. V25–V32. External Links: Document Cited by: §1, §1, §3.3. Romera-Paredes et al. (2024) B. Romera-Paredes, M. Barekatain, A. Novikov, M. Balog, M. P. Kumar, E. Dupont, F. J. R. Ruiz, J. S. Ellenberg, P. Wang, O. Fawzi, P. Kohli, and A. Fawzi Mathematical discoveries from program search with large language models. Nature 625 (7995), p. 468–475. External Links: Document Cited by: §1. Sheng et al. (2025) H. Sheng, X. Wu, X. Si, J. Li, S. Zhang, and X. Duan Seismic foundation model: a next generation deep-learning model in geophysics. Geophysics 90 (2), p. IM59–IM79. External Links: Document Cited by: §1. Spitz (1991) S. Spitz Seismic trace interpolation in the F-X domain. Geophysics 56 (6), p. 785–794. External Links: Document Cited by: §1, §1. Sun et al. (2026) B. Sun, W. Guo, Z. Yu, and L. Yang Self-evolving scientific agent discovers generalizable physically-reasoned fluid control. arXiv preprint arXiv:2606.08405. Cited by: §1. Trad (2009) D. Trad Five-dimensional interpolation: recovering from acquisition constraints. Geophysics 74 (6), p. V123–V132. Cited by: §1. Wang et al. (2014) B. Wang, R. Wu, Y. Geng, and X. Chen Dreamlet-based interpolation using POCS method. Journal of Applied Geophysics 109, p. 256–265. Cited by: §1, §2.1. Wang et al. (2019) B. Wang, N. Zhang, W. Lu, and J. Wang Deep-learning-based seismic data interpolation: a preliminary result. Geophysics 84 (1), p. V11–V20. External Links: Document Cited by: §1. Wu and Bai (2018) J. Wu and M. Bai Adaptive rank-reduction method for seismic data reconstruction. Journal of Geophysics and Engineering 15 (4), p. 1688–1703. External Links: Document Cited by: §1, §3.3. Xu et al. (2026) Y. Xu, W. Wang, S. Yu, and J. Ma Unsupervised anti-aliasing interpolation of regularly sampled seismic data via soft dynamic time warping divergence distance. Petroleum Science. Note: In press External Links: Document Cited by: §1. Xu et al. (2024) Y. Xu, S. Yu, L. Dong, and J. Ma Dealiased seismic data interpolation by dynamic matching. Geophysics 89 (5), p. V361–V376. External Links: Document Cited by: §1. Yang et al. (2026) M. Yang, L. Guo, and L. Yang Harness in-context operator learning with chain of operators. arXiv preprint arXiv:2606.12318. Cited by: §1. Ye et al. (2024) H. Ye, J. Wang, Z. Cao, F. Berto, C. Hua, H. Kim, J. Park, and G. Song ReEvo: large language models as hyper-heuristics with reflective evolution. In Advances in Neural Information Processing Systems, Vol. 37, p. 43571–43608. External Links: Link Cited by: §1. Yu et al. (2019) S. Yu, J. Ma, and W. Wang Deep learning for denoising. Geophysics 84 (6), p. V333–V350. External Links: Document Cited by: §1. Yu and Yang (2026a) Z. Yu and L. Yang Agentic symbolic search: characterizing pdes beyond hand-crafted expressions, meshes, and neural networks. arXiv preprint arXiv:2606.20467. Cited by: §1. Yu and Yang (2026b) Z. Yu and L. Yang Evolutionary ensemble of agents. arXiv preprint arXiv:2605.09018. External Links: 2605.09018, Document, Link Cited by: §1. Zhang et al. (2016) D. Zhang, Y. Chen, W. Huang, and S. Gan Multi-step damped multichannel singular spectrum analysis for simultaneous reconstruction and denoising of 3D seismic data. Journal of Geophysics and Engineering 13 (5), p. 704–721. Cited by: §1, §2.1. Zhang et al. (2015) H. Zhang, X. Chen, and H. Li 3D seismic data reconstruction based on complex-valued curvelet transform in frequency domain. Journal of Applied Geophysics 113, p. 64–73. Cited by: §1, §2.1.