# Optimal experimental design for full waveform inversion using a wavenumber sampling criterion-part 1: 2-D to methodological development

Amine Abdellaziz, Romain Brossier, Ludovic Métivier, Édouard Oudet

► To cite this version:

Amine Abdellaziz, Romain Brossier, Ludovic Métivier, Édouard Oudet. Optimal experimental design for full waveform inversion using a wavenumber sampling criterion-part 1: 2-D to methodological development. Geophysical Journal International, 2025, 240 (3), pp.1429-1459. 10.1093/gji/ggae454. hal-05336522

HAL Id: hal-05336522

https://hal.science/hal-05336522v1

Submitted on 29 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal experimental design for full waveform inversion using a wavenumber sampling criterion—part 1 : 2-D to methodological development

Amine Abdellaziz<sup>*</sup>,<sup>1,2</sup> Romain Brossier,<sup>1</sup> Ludovic Métivier<sup>1,2</sup> and Édouard Oudet<sup>2</sup>

Accepted 2024 December 18. Received 2024 December 12; in original form 2024 October 4

## SUMMARY

Due to the complexity of seismic surveys, different strategies have been investigated to optimize characteristics of the acquisition in order to fulfil various criteria, most intuitively those related to the operational and economical aspects of the survey. However the question can also be approached from an imaging perspective. Andrzej Kikjo's seminal work on the optimal design of seismic stations for earthquake localization, dating back to the seventies, has encouraged several works related to optimal acquisition in exploration geophysics. A common choice for the optimality criterion is the conditioning of the Hessian operator. In our work, we have chosen our optimality criterion to be the regularity of the wavenumber sampling at a specified target: using the theory of diffraction tomography we can estimate the wavenumber sampling at a target point and we note that for a regularly spaced acquisition the corresponding sampling is not regular. In our target-oriented approach, we want to find the best positions of sources and receivers to improve the regularity of this sampling. We develop our work in the context of full waveform inversion, although our approach is more general. Using a simple formula for the estimation of the wavenumber points we are able to get a geometric characterization of the accessible wavenumber space that helps us define its envelope. We describe the problem as finding the geometry for which the sampling points are regular in terms of spacing inside this envelope and to that end we incorporate in our workflow the notion of centroidal Voronoi tessellation, a tool from computational geometry that links the regularity of the sampling to the uniformity of the domain tessellation. To assess the uniformity of the tessellation we utilize an energy function the minimal values of which provide regular samplings. By combining this energy function with a map that estimates the wavenumber content of an acquisition we define a new objective function that assess the quality of such acquisition with regards to our regularity criterion. We then mathematically express our problematic as a minimization problem which we solve using local optimization strategies. In this article we focus on the theoretical aspect of our method and validate it through numerical examples of varying complexities in the 2-D cases: when comparing optimized acquisitions with regular ones the results show a consistent local improvement of the quality of the inversion with the optimized geometries at the target and motivates the exploitation of the method for target-oriented applications such as seismic monitoring. The subsequent extension to 3-D and application to more realistic industrial examples will be discussed in a companion paper.

**Key words:** Inverse theory; Controlled source seismology; Crustal imaging; Seismic instruments; Waveform inversion.

## 1 INTRODUCTION

Seismology is used to study the structure of the subsurface. To get information about the latter we generally cannot rely on *in situ* measurements.

For example drilling wells is not conceivable for large scale surveys since it only provides local information on the subsurface, or for archaeology where it would damage the target. Moreover drilling is expensive and despite advances in

<footnote n="1"><sup>1</sup>Univ. Grenoble Alpes, ISTERRE, F-38000 Grenoble, France. E-mail: mohamed-amine.abdellaziz@univ-grenoble-alpes.fr</footnote>

<footnote n="2"><sup>2</sup>Univ. Grenoble Alpes, CNRS, LJK, F-38000 Grenoble, France</footnote>
technology the depth of exploration is limited (the deepest drilling operation managed to go slightly beyond 12 km; Smythe et al. 1994). Therefore we need to rely on indirect physical measurements. Such measurements include gravity, electromagnetic (EM) waves and mechanical waves. Gravity survey highlights density anomalies that can be related to geological features and is a low-resolution imaging technique (Blakely 1995). The variation in structure also affects the conductivity of the subsurface so EM waves are used to detect them. However, high-frequency propagative EM waves have a very small penetration depth of few tens of metres maximum. At low frequency, diffusive EM wave can reach deeper depths down to few (tens) of kilometres, but with very low resolution. In a similar manner, by measuring mechanical waves we can infer the parameters of the subsurface (velocity, density, attenuation and anisotropy) that affect the propagation. Those waves have a deeper penetration and therefore are used at all scales depending on their frequency-content (global, regional, crustal, near-surface).

To characterize the structure and mechanical parameters of a region inside the Earth the idea is either to rely on natural seismic sources (earthquake, volcanic activity,...) or to position devices that we call active sources that generate perturbations and emit mechanical waves propagating through the medium. Those waves are recorded by other devices called receivers. The recorded time-series of those vibrations are assembled into seismograms. From the seismograms, we want to find properties of the medium that best describe them. This is a nonlinear inverse problem which is ill-posed: the solution is highly sensitive to changes in the data and the choice of the initial guess. Given the central role that sources and receivers play in this process it is natural to ask the question: what is the effect of the positions on the reconstruction? Can we determine criteria to select the best positions for our experiment?

Optimal experimental design (OED) is a field of study that originally comes from statistics. It provides a framework within which to obtain high-quality designs by optimizing a number of criteria generally related to the variance of the sought-after parameters or other statistically relevant properties (Smucker et al. 2018). First theoretically developed by Smith (1918), it has been used in many fields where experiments are conducted: examples include clinical trials for the development of medical drugs (Khinkis et al. 2009), the estimation of heat and moisture for material science (Berger et al. 2017) or optimal control (Pronzato 2008). Optimal design has first been used in geophysics by Kijko (1977) for the localization of earthquakes. Using Geiger's method (see Karasözen & Karasözen 2020, for a review), a linearized earthquake localization method, and assuming a knowledge of the subsurface model, the method developed by Kijko aims at finding the positions for seismic stations that ensure the best localization of all earthquake foci and origin times. In his work, he looks for designs which minimize the area of the ellipse volume of the location error generated from seismic random errors. The criterion is also called D-optimal criterion and corresponds to maximizing the determinant of the information matrix (an approximation of the Hessian) which has the effect of improving its conditioning. In OED near-optimal designs are considered adequate because they still provide an improvement. Therefore solutions are computed using numerical heuristics that are either specifically developed for OED, for example the DET-MAX algorithm (Rabinowitz & Steinberg 1990; Mitchell 2000), or more general algorithms such as simulated annealing (Hardt & Scherbaum 1994).

In passive seismology we do not have control over the sources since we rely on phenomena such as earthquakes to record information. Thus, we can only control the positioning of receivers. In the context of geophysical exploration, we use controlled seismic sources and thus we decide where to position them, and we have more degrees of freedom to design the acquisition. Besides, depending on the type of application, we have access to different choices of equipment to use in the survey. This in turn affects the geometry of acquisition.

For onshore applications the type of sources depends on the scale of the survey. For near-surface characterization, a sledgehammer hitting a metal plate on the ground produces sufficient energy to illuminate the target zone. For exploration scale, more powerful sources are needed for the wave to penetrate deeper. Seismic vibrators mounted on trucks are often used, and sometimes explosive sources are used. For the receivers, geophones recording three components of the particle field velocity are used. Formerly the receivers had to be connected with wires which limited the possibilities of deployment; the receivers were placed in lines and this restricted the acquisition to grid-like geometries. Nowadays the technology moved to wireless nodes which lifts such restrictions. Moreover lightweight nodes have been developed (Ourabah et al. 2020) that can easily be transported and the use of drones have been proposed to cover usually inaccessible areas such as foothills (Pagliccia et al. 2018).

For offshore applications the most common sources are airguns, devices that expel into the water high-pressure air stored in a chamber. In water, hydrophones are used to measure pressure and they are placed inside cables towed by streamer boats. Sources and receivers can be towed by the same boat or distinct ones. In both cases the geometry of the acquisition is limited by the trajectory taken by the boat. Ocean bottom node (OBN) is a technology located on the seafloor which has gained traction as it consists of a three components geophones and a hydrophone and it allows access to larger offsets. However they are operationally costly since they need to be deployed before the acquisition and retrieved later; those operations represent a significant amount of survey time. To alleviate those limitations new technologies are developing where OBNs are mounted in autonomous underwater vehicles (Mancini et al. 2023) and surface floating drones (Amelin & Yakovlev 2024) to automate their deployment and retrieval. Like onshore nodes, OBNs facilitate the flexibility of the acquisition in the sense that we can design arbitrary geometries where the positions of sources and receivers are independent from each others.

The dimensions of the survey vary depending on the application. The largest campaigns can cover areas of hundreds of square kilometres, deploy tens of thousands of receivers and last for several months. It thus involves large crews and considerable amount of equipment. Moreover seismic surveys face various constraints like environmental requirements dictated by local authorities (for instance restricted time window for shooting sources in order to protect preserved natural zones) or deployment limits (for instance avoid obstacles such as wind farms or maritime routes). An example of large and complex survey campaign can be found in Baudo et al. (2023).

A seismic survey is therefore operationally heavy and costly. Given the scale and complexity of seismic surveys, it is only natural that efforts are done to find acquisitions adequate for specific purposes.

A factor of interest in seismic acquisition, especially for operators, is the cost of the survey. This cost is ultimately controlled by the
time of the acquisition and falls under operational considerations. Morrice *et al.* (2001) proposes a mathematical model to optimize the design of land acquisition to minimize operational costs. He considers an orthogonal acquisition (the source and receivers lines are orthogonal) and wants to find the optimal design parameters (for instance source and receiver stations spacing, line spacings) that minimizes an objective function evaluating the cost of the survey (including costs of the crew, of the survey permit, drilling, equipment and operations). He also adds geophysical constraints (requirements in terms of bin size, fold, offset and azimuth) and operational constraints (limit of the equipment available, maximum number of shots acquired by day). He expresses the problem as a mathematical programming model. Without such formal description of optimal design, various heuristic strategies are employed in the industry to reduce the operational costs of a survey. We can cite the choice of spiral trajectory for the seismic source boat in marine survey (Ridyard *et al.* 2023) to reduce unnecessary turns and manoeuvres and thus saving time. To also reduce the time of the acquisition we can decide to fire several sources simultaneously to bring down the length of the survey campaign and later use deblending to separate the origins of the arrivals (Wang *et al.* 2024). Another strategy which is proposed to speed up acquisition and decrease equipment use is simply the decimation of the acquisition (Pagliccia *et al.* 2018). For that the application of compressive sensing in seismic imaging is explored, which relies on signal processing and enables sampling at sub-Nyquist rates without loss of information (Baraniuk & Steeghs 2017).

While the cost of acquisition survey is heavy, the acquisition remains crucial for imaging. Therefore another approach for the optimal design consists in considering geophysical criteria instead of operational ones. Liner *et al.* (1999) proposes to find land acquisition parameters (source and receiver stations spacing, source and receiver lines spacing, inline and cross-line offsets) to achieve desirable characteristics (target values for bin-size, offset and fold). Here we are looking for acquisition designs that fit geophysical requirements. So if we know which geophysical characteristics are desirable in the design then as long as those requirements are quantifiable and maximizable then we can optimize the acquisition with regards to those criteria. Stork (2019, 2021) explores the optimization of sources and receivers locations to improve the signal-to-noise ratio by avoiding the most noisy areas. A promising alternative to conventional acquisitions is distributed acoustic sensing (DAS), and Fichtner & Hofstede (2022) propose a simple algorithm for optimizing DAS acquisition designs. For that, they randomly generate different cable layouts and keep the configuration that best fit a certain criterion, such as the homogeneity of the offset coverage across azimuth bins. The method is highly adaptable and can be combined with various optimization criteria. Finally another approach chooses, similarly to the earthquake localization problem, the criterion for optimization to be the conditioning of the imaging operator associated to the inverse problem. A benefit measure is the eigenvalue spectrum of the approximate Hessian (Curtis & Snieder 1997). Such a criterion is computationally expensive and in the literature approximate measures have been explored relying on the trace of the matrix (Curtis 1999) or its determinant (Krampe *et al.* 2020). This connects to the conditioning of the operator that the D-optimality tries to improve.

We have seen so far that we can address the optimal design of acquisition from two perspectives: either from the point of view of the survey operators (who consider economical and operational constraints) or from the point of view of the imaging (knowing how the data will be used afterwards). In our work we focus on

the optimal design for the imaging of a given target based on a full
waveform inversion (FWI) approach, which is an imaging method
relying on the iterative minimization of the misfit between observed
and calculated data. It has been introduced in the 80s and is now
routinely applied both in academy and industry (Virieux & Operto
2009; Virieux *et al.* 2017).

We propose an alternative approach for the optimal design of ac-
quisition geometry for exploration geophysics. We consider fixed-
spread surface acquisitions (each source is visible to every receiver,
so there is no displacement of the receivers between shots) where
the positions of sources and receivers are independent from each
other (arbitrary geometry). We fix a target in the subsurface and we
optimize the positions of the sources and receivers in order to im-
prove the imaging of said target. To that end the selected criterion is
the wavenumber sampling provided by the acquisition at the target:
we want to find the acquisition geometry for which this sampling
is as regular as possible. This approach helps to mitigate problems
related to undersampling and aliasing, which can degrade the qual-
ity of seismic images and disrupt subsurface imaging (Biondi 2006,
chapter 8).

Our new method relies on three complementary tools :

(1) A function that maps the acquisition to the wavenumber sam-
pling at the target. This sampling is estimated using the diffraction
tomography approximation.

(2) A Voronoi tessellation tool from computational geometry that provides a way to quantify the regularity of the sampling.

(3) An objective function that we define using the two previous tools in order to assess the quality of the acquisition. We express our optimal design problem as a minimization problem and solve it using methods of mathematical optimization.

We will expose our work in two papers: in this first paper we will lay the theoretical foundation of our approach and illustrate it in the context of 2-D examples. In a second paper, we will present the extension of this method to 3-D application, and insert it in more realistic survey design scenarios by taking into account deployment constraints imposed on the acquisition.

## 2 METHODOLOGY

In this section we present our method. At first, we present how the wavenumber content associated with an acquisition is estimated given a target. We will study its effect on the imaging and characterize its distribution for surface acquisition. Then we discuss how to take into account the effect of the heterogeneity of the medium in the illumination. Afterwards we introduce our Voronoi tessellation tool and show how to regularize the sampling by regularizing the tessellation. Finally by combining the Voronoi tessellation and the mapping from acquisition to the wavenumber space we develop an approach for the optimization of the acquisition with regards to the regularity sampling in the wavenumber space.

### 2.1 Acquisition to wavenumber space mapping

In a 2-D space let **x** be the position of the target to image and consider an acquisition of $N_s$ sources and $N_r$ receivers. The FWI gradient $\nabla f(\mathbf{x})$ can be computed as the cross-correlation between incident wavefields $u_s(\mathbf{x})$ coming from the sources and adjoint wavefields $\lambda_r(\mathbf{x})$ coming from the receivers (Virieux & Operto 2009). From diffraction tomography (Wu & Toksöz 1987), assuming plane wave propagation from the source $u_s(\mathbf{x}) \approx \exp(i\mathbf{p}_s \cdot \mathbf{x})$ and the receiver
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic representation of eq. (2). Source is denoted by **s** and receiver by **r**. Rays connecting the devices to the target **x** are considered as straight lines. The directions of those rays are given by unit vectors **p**<sub>s</sub> and **p**<sub>r</sub>, the sum of which yields the wavenumber **k**.</figcaption></figure>

$\lambda_r(\mathbf{x}) \approx \exp(i\mathbf{p}_r \cdot \mathbf{x})$, it can be shown that the cross-correlation becomes a sum over the sources and the receivers of products of plane waves and the gradient can be expressed as an inverse Fourier summation (Sirgue & Pratt 2004):

$$ \nabla f(\mathbf{x}) \approx \sum_{s \in S} \sum_{r \in R} \exp(ik_0 (\mathbf{p}_s + \mathbf{p}_r) \cdot \mathbf{x}). \quad (1) $$

The wavenumber content of the approximated gradient is generated from each source-receiver couple and is linked to the geometry of the acquisition through the eq. (1), from which we can write the resulting wavenumber vector

$$ \bar{k}(\mathbf{s}, \mathbf{r}) = k_0 (\mathbf{p}_s + \mathbf{p}_r), \quad (2) $$

where **p**<sub>s</sub> and **p**<sub>r</sub> are the unit direction vectors of the rays connecting the source **s** and the receiver **r** to the target **x** (Fig. 1). The reference wavenumber of the two plane waves is given by $k_0 = 2\pi f_0/c_0$, where $f_0$ is the reference frequency and $c_0$ is the velocity at the target.

All along our study, we make the assumption that the acquisition is fixed-spread. This means that each source sees each receiver. In practice from eq. (2), we get a sampling of $N_s \times N_r$ wavenumbers. However from an optimal design point of view, the number of degrees of freedom is still the total number of sources and receivers, namely $N_s + N_r$. We define a mapping from the acquisition to the wavenumber space

$$ G : \begin{array}{ccc} D_s^{N_s} \times D_r^{N_r} & \rightarrow & W^{N_s \times N_r} \\ (S, R) & \mapsto & K = G(S, R). \end{array} \quad (3) $$

The map **G** takes as input a set of sources $S = (\mathbf{s}_i)_{i=1,...,N_s}$ and receivers $R = (\mathbf{r}_j)_{j=1,...,N_r}$ and computes the set $K = (k(\mathbf{s}, \mathbf{r}))_{\mathbf{s} \in S, \mathbf{r} \in R}$ of associated wavenumbers using eq. (2). The sets $D_s$ and $D_r$ define the domains of possible positions of, respectively, the sources and receivers. The wavenumber domain is represented by the subset $W \subset \mathbb{R}^2$ of all possible samples of the wavenumber space. Because $\mathbf{p}_s$ and $\mathbf{p}_r$ are unit vectors in (2), we deduce that the wavenumbers will not exceed a magnitude of $2k_0$, therefore the domain $W$ will be included inside a closed disc $\bar{D}(O, 2k_0)$ in the wavenumber space of centre the origin $O = (0, 0)$ and radius $2k_0$.

To illustrate our aim we conduct a small imaging experiment. In Fig. 2 we have a $V_p$ velocity model of dimensions 1 km × 1 km and background velocity of 2000 m s<sup>-1</sup>. At the centre of the model is a small Gaussian anomaly of approximately 50 m of diameter with a peak velocity at 2500 m s<sup>-1</sup>. By defining this anomaly as the target, we want to study the effect of the geometry, alternatively the wavenumber sampling in the inverse Fourier transform of the

<figure><img src="image_3.png" /><figcaption>**Figure 2.** True model of the experiment in Fig. 3. The model is of dimensions 1 km × 1 km and has a homogeneous background of 2000 m s<sup>-1</sup> with a Gaussian anomaly at the centre. The anomaly represents the target we want to image using FWI and taking the homogeneous background as initial model.</figcaption></figure>

approximation in eq. (1), on the reconstruction of the FWI gradient. For that we set the model with anomaly as the true model and the initial one is the homogeneous background. To test the effect of the acquisition, we consider different geometries. All acquisitions have 32 sources and 32 receivers. We define three different settings: (a) a circle around the target, (b) a cross-hole acquisition and (c) a surface acquisition. In Fig. 3 we show for each acquisition the geometry of the acquisition (left panel) along with the wavenumber sampling (centre panel) and the normalized FWI gradient (right panel).

In the first row (a) the acquisition is a perfect circle around the target. The sources and receivers are alternating and regularly spaced. Because of the circular layout we have a full angular coverage at the target and therefore the wavenumbers sample the whole area of the disc $\bar{D}(O, 2k_0)$ in the wavenumber space. Thanks to the dense sampling (1024 = 32 × 32 sample points) and the complete angle coverage, the perturbation model provided by the FWI gradient is well focused at the position of the anomaly in the true model. It corresponds to a low-pass filtered version of the target due to the band-limited frequency content of the data and would correspond to a good reconstruction of the image.

In the second case, Fig. 3 row (b), the acquisition mimics a cross-hole setting with sources and receivers positioned in different wells on each side of the target. They illuminate the target from two different sides, so the $\mathbf{p}_s$ and $\mathbf{p}_r$ vectors have opposite horizontal directions. This has a cancelling effect on the $k_x$-component of the wavenumbers so the sampling consists of points with small $k_x$ values. When comparing the gradient with the circular case (a) we see that we have a comparably good reconstruction of the vertical profile, contrary to the horizontal reconstruction.

In the third row (c) we consider a surface acquisition with alternating and regularly spaced sources and receivers. Here the target, located below the acquisition, is illuminated only from above so the
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Acquisitions with three different settings along with their wavenumber content at the target and L2-normalized first FWI gradient. The acquisitions consider 32 sources and 32 receivers which yields a sample of 1024 points in the wavenumber space. Here we have set the reference wavenumber $k_0 = 2$. This simple experiment highlights the impact of the acquisition geometry on the wavenumber sampling and FWI imaging. The FWI gradient has been computed using the TOYxDAC_TIME program (Yang *et al.* 2018).</figcaption></figure>

$p_s$ and $p_r$ vectors will point downwards and generate wavenumbers with positive $k_z$-component only. This limited illumination also restricts the sampling of the vertical wavenumber space to a positive range. This creates a band-pass filtered wavenumber reconstruction vertically with the resulting Cardinal-Sine-like shape on the vertical section. On the other hand, the horizontal update shows a low-pass filtered version of the anomaly, thanks to the symmetric horizontal coverage in $k_x$.

From these three canonical examples, we can see that the acquisition geometry dramatically affects the reconstruction quality, with a direct link with the wavenumber sampling. Preferably, the covered area of the wavenumber space would need to be as large as possible and the distribution of wavenumbers inside this area as regular in terms of spacing as possible.

In our work, we focus on the surface acquisition case. A first step to the optimization of the coverage is the study of the wavenumber
<figure><img src="image_5.png" /><figcaption>Highlighted circular arcs</figcaption></figure>

**Figure 4.** Illustration of the circular arc characterization of the wavenumber domain $W$. We consider a surface acquisition (in the left panel) with dense receivers and few sources. We note that the wavenumber cloud (in the right panel) is composed of circular arcs. Each circular arc and its centre (crosses) are of the same colour as the source from which it has been generated. Note that the leftmost (resp. rightmost) arc corresponds to the rightmost (resp. leftmost) source in the acquisition and the centres of arcs form also a circular arc pattern.

cloud and more precisely the domain $W$ of all possible wavenumbers. We can describe the wavenumber cloud as a union of circular arcs. Indeed, we want to find the trajectory made by all the wavenumbers associated with a single source. For that we fix the position of the source $\mathbf{s}$. Then the vector $\mathbf{p}_s$ will be fixed and only $\mathbf{p}_r$ will vary. Because both are unit vectors then from eq. (2) the wavenumbers associated to the source will be a circular arc. By doing this analysis for every source in the acquisition, the cloud of points will be a union of arcs.

Mathematically we study the wavenumbers using an angular characterization and therefore take advantage of trigonometric identities. As shown in Fig. 1 the direction vectors $\mathbf{p}_s$ and $\mathbf{p}_r$ can be expressed using the illumination angles $\phi_s$ and $\phi_r$ they make with the downward pointing vertical axis of direction $\bar{n}$. The coordinates of those vectors are given by $\mathbf{p}_s = [\sin \phi_s, \cos \phi_s]^T$ and $\mathbf{p}_r = [\sin \phi_r, \cos \phi_r]^T$. We can then get an angular expression of formula (2)

$$k(\mathbf{s}, \mathbf{r}) = k_0 \left[ \frac{\sin \phi_s + \sin \phi_r}{\cos \phi_s + \cos \phi_r} \right]. \qquad (4)$$

Let us again fix the position of the source, giving a fixed angle $\bar{\phi}_s$ and let the location of the receiver vary which results in the other angle $\phi_r$ varying between $\phi_{min}$ and $\phi_{max}$. Here $[\phi_{min}, \phi_{max}]$ denotes the illumination range the directional vectors $\mathbf{p}_s$ or $\mathbf{p}_r$ can make with the vertical and is related to the right and left limits of the acquisition. From eq. (4) we would have for any $\phi_r \in [\phi_{min}, \phi_{max}]$

$$k(\mathbf{s}, \mathbf{r}) = k_0 \mathbf{p}_s + \left[ \frac{k_0 \sin(\phi_r)}{k_0 \cos(\phi_r)} \right]. \qquad (5)$$

Consequently, the wavenumbers associated with the source $\mathbf{s}$ create a circular arc of radius $r = k_0$, centre $c = k_0 \mathbf{p}_s = k_0 (\sin \bar{\phi}_s, \cos \bar{\phi}_s)$ and length $l = k_0 (\phi_{max} - \phi_{min})$. The same analysis can be done for a receiver $\mathbf{r}$, with its fixed angle $\bar{\phi}_r$ replacing $\bar{\phi}_s$ and letting this later vary. Considering a deployment range on the surface that produces minimum and maximum illumination angles $\phi_{min}$ and $\phi_{max}$, the

domain $W$ can be defined as

$$ \begin{aligned} W &= \bigcup_{\phi_s = \phi_{min}}^{\phi_{max}} \left\{ k_0 \mathbf{p}_s + k_0 \left[ \frac{\sin(\phi_r)}{\cos(\phi_r)} \right] : \phi_{min} \le \phi_r \le \phi_{max} \right\} \\ &= \bigcup_{\phi_s = \phi_{min}}^{\phi_{max}} \mathrm{Arc}_{\phi_s}, \end{aligned} \qquad (6) $$

where $\mathrm{Arc}_{\phi_s}$ denotes the circular arc of radius $r = k_0$ and centre $c = k_0 (\sin \phi_s, \cos \phi_s) = k_0 \mathbf{p}_s$, generated by angles $\phi_r$ spanning the illumination range $[\phi_{min}, \phi_{max}]$.

The domain of sampled wavenumbers $W$ has therefore an envelope constituted by the union of three arcs: the two side circular arcs $\mathrm{Arc}_{\phi_{min}}$ and $\mathrm{Arc}_{\phi_{max}}$, in addition to the bottom arc of centre the origin, radius $r = 2k_0$ and angle $\phi_s \in [\phi_{min}, \phi_{max}]$. The bottom arc is composed of wavenumbers with the largest amplitude which are associated with the sources and receivers located at the same position:

$$k(\mathbf{s}, \mathbf{r}) = 2k_0 \left[ \frac{\sin \phi_s}{\cos \phi_s} \right] \quad \text{for } \phi_s \in [\phi_{min}, \phi_{max}]. \qquad (7)$$

Fig. 4 illustrates this characterization. We have on the left panel a surface acquisition with a dense receiver line (50 devices) and a small number of sources (only 8). On the right panel, we show the wavenumber sampling at the target. We can clearly see that the wavenumber cloud is constituted of eight arcs, each one associated with a source in the acquisition. To distinguish between them, a distinct colour is assigned to each source and its corresponding arc. In the wavenumber space, we also plot the centres of arcs and we note that they form a circular trajectory. This can easily be seen when we take a look at the set $\{k_0 \mathbf{p}_s = k_0 [\sin \phi_s, \cos \phi_s]^T | \phi_s \in [\phi_{min}, \phi_{max}]\}$ of the possible centres of arcs which itself defines an arc.<fnref n="1" />

From this geometric characterization, we can better understand the influence of the frequency $f_0$ and the offset on the sampling.

In Fig. 5 we study the effect of the frequency on the wavenumber coverage. On the left panel, we have a surface acquisition with

<footnote n="1">¹Please note that these centres are presented here for the sake of explaining the wavenumber coverage. They do not belong to the wavenumber cloud associated with the acquisition.</footnote>
Effect of the frequency

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Effects of the frequency on the wavenumber sampling. On the surface we have a regularly spaced acquisition (left panel) and we generate the wavenumbers through the map $G$ using three different frequencies $f_0$ (right panel). Increasing the frequency has the effect of increasing the radius of the circular arcs, therefore the size of the sampled area.</figcaption></figure>

Effect of the offset

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Effects of the offset on the wavenumber sampling. On the surface we deploy two acquisitions (left panel) with different offsets and we generate the corresponding wavenumbers (right panel) through the map $G$ along with the envelope of their respective domain $\mathcal{W}$. Insolid line we also plot the envelope of the largest area which we can theoretically cover in a homogeneous medium if we had access to an infinite deployment range on the surface. This represents an upper bound on the covered area by the wavenumber sampling.</figcaption></figure>

a target in the subsurface. The medium is homogeneous and on the right panel we compute the wavenumber cloud using the map $G$ for different frequencies $f_0$ (namely 0.5 Hz, 1 Hz and 1.5 Hz). We note that the area of the point cloud is correlated with the frequency: a larger frequency yields larger coverage. This is easy to grasp since increasing the frequency increases the value of the reference wavenumber $k_0$ and consequently the magnitude of the wavenumbers. From the geometric characterization provided by (6), increasing $k_0$ has the effect of increasing the radius $r = k_0$ of the circular arcs, thus the size of the domain $\mathcal{W}$.

In Fig. 6 we study the effect of the offset on the wavenumber sampling. The left panel shows two acquisition geometries with respectively large offsets (in blue) and small offsets (in orange, half the maximum offset of the blue one). On the right panel, we plot the wavenumber clouds and their envelopes. We can see that the offset also affects the size of $\mathcal{W}$ but in a different way than the frequency: we have here a lateral expansion of the size. This is due to the link between the offset and the maximum/minimum illumination angles $\phi_{\max}$ and $\phi_{\min}$. When changing the offset range, the length of the circular arcs are affected: a larger offset range is associated with a larger $(\phi_{\max} - \phi_{\min})$ spectrum and thus larger arcs. In a homogeneous medium, the illumination angle is bounded by $\pi/2$, so the size of the wavenumber envelope is also bounded. For the

sake of illustration, we plot in red the upper bound of the envelope corresponding to an acquisition with an infinite offset length that would allow to reach an illumination angle $\phi_{\max} = \pi/2$ from both sides of the acquisition.

In our methodology, we consider that we fix the frequency $f_0$ and the maximum offset allowed by the acquisition. Therefore, the envelope area in the wavenumber space is fixed by those two settings. We also fix the number of available $N_s$ sources and $N_r$ receivers, fixing the number of sampling points in the wavenumber domain. Our optimization will then influence the sole remaining free parameters: the positions of sources and receivers. In Fig. 7 we show how different surface geometries with the same number of sources and receivers get different wavenumber sampling. On top a regularly spaced surface acquisition is considered. At the target point, we also show the illumination ($\mathbf{p}_s$ and $\mathbf{p}_r$ vectors). We can note a duality between the geometry and the illumination: regularity in the geometry (positions) does not translate into a regular illumination (angles). On the right panel, we show the wavenumber sampling along with the envelope of the domain $\mathcal{W}$. The irregular sampling in angles yields a wavenumber sampling which is also irregular inside the envelope: the points tend to concentrate at the tips of the envelope whereas the centre is sparsely sampled.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Comparison between a regular and optimized surface acquisition in terms of wavenumber sampling. In the acquisitions (on the left panel) we have 10 sources and 11 receivers positioned along a 4 km line and a target at 1 km depth. At the target we plot the illumination ($\mathbf{p}_s$ and $\mathbf{p}_r$ vectors) and in the wavenumber space (right panel) the sampling. Because the two acquisitions share the same deployment range and frequency $f_0$, the domain $\mathcal{W}$ is the same for both. Also the samplings have the same number of points because the two acquisitions have the same number of sources and receivers.</figcaption></figure>

Our objective is therefore to find an acquisition geometry for which the wavenumber sampling is more regularly distributed. In Fig. 7, we show in the bottom row an example of optimized geometry and the associated wavenumber plot. We can clearly see that for this particular geometry, the illumination and sampling are comparatively regular, at the expense of an irregular positioning at the surface. This represents the kind of results we would like to achieve.

This 'optimal' geometry has been obtained with our method, which we will present in the next part of this section. Before discussing the inner working of our algorithm, we will first discuss how to consider the heterogeneity of the medium in the estimation of the wavenumbers.

## 2.2 Heterogeneous medium

In a target-oriented approach we need to take into account as much as possible the prior knowledge of the medium properties we have access to. Although we have described the travelling rays using straight-lines (corresponding to a plane-wave approximation in a homogeneous medium) to evaluate the wavenumber content, our approach is valid in any medium for which the direction vectors $\mathbf{p}_s$ and $\mathbf{p}_r$ can be calculated.

In order to deal with heterogeneous media, we rely on the ray theory as an asymptotic high frequency approximation to estimate the direction vectors. The traveltime $T(x)$ of a wave originating from a point source in the medium obeys an eikonal

equation

$$ \|\nabla T(x)\| = \frac{1}{v(x)}, \qquad (8) $$

with $v(x)$ representing the velocity of the medium at the point $x$. In isotropic media the rays are parallel to the slowness vector given by the gradient $\nabla T(x)$ of the traveltime (Le Bouteiller 2018). Hence the direction of the ray connecting the source (resp. receiver) to the target is co-linear to the gradient of the arrival times $T(x)$ of the wave field originating from the source (resp. receiver) and evaluated at the target. This means that by computing the values of the gradient $\nabla T(x)$ at the target $x$ we can get the directional vector $\mathbf{p}_s$ (resp. $\mathbf{p}_r$) through equation

$$ \mathbf{p}_s = \frac{\nabla T(x)}{\|\nabla T(x)\|}. \qquad (9) $$

Fig. 8 illustrates this process.

To illustrate the contribution of the heterogeneity to the wavenumber sampling, we show in Fig. 9 an example of acquisition in a heterogeneous medium along with its wavenumber content. The acquisition is composed of 10 sources and 11 receivers spanning a deployment area of 8 km on the surface. The target is situated at the centre of the model at a depth of 1 km. The medium has a 1-D velocity model ranging from 1500 m s⁻¹ at the surface to 4000 m s⁻¹ at 2 km depth. Due to the large offset, we can see that we have access to diving waves and this means that we are able to widen the illumination. As a result on the right panel we can see that we have a wavenumber sampling with a significantly larger coverage than if we only had a homogeneous model. Despite the particular shape of
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Illustration of ray-tracing. In this example, we have a heterogeneous velocity model (top colour map) and would like to find the ray connecting the source (triangle) to the target (cross). For that we compute the arrival times of the wave propagating from the source (bottom colour map) by numerically solving the eikonal eq. (8) and we compute the gradient at points of the model (arrows). The direction $p_s$ follows the same direction as the gradient at the target. The arrival times and the gradient have been numerically computed using the HFM program (Mirebeau & Portegies 2018).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Example of wavenumber sampling in a heterogeneous medium. We have a surface acquisition (left panel) on a deployment range of 8 km and the velocity model is 1-D. The target is at depth 1 km. Also shown in the panel are the rays connecting the sources and receivers to the target and the illumination represented by direction vectors $p_s$ and $p_r$. We plot the corresponding wavenumber cloud (right panel) through the map $G$. Note the particular shape of the domain $W$ compared with the envelope in a homogeneous medium (Fig. 7). Despite the difference in shape, both domains are characterized using circular arcs through eq. (6).</figcaption></figure>

the wavenumber envelope, the circular arc characterization of the wavenumber cloud is still valid: the arcs are just larger and thus they intersect. Indeed with the undershooting we have the illumination range [$\phi_{\min}, \phi_{\max}$] which exceeds the upper bound range [−π/2, π/2] accessible when considering straight rays. In turn, the length of the

arcs exceed $k_0\pi$ and consequently the two side arcs of the envelope intersect at the origin.

As it will be discussed later, it is advantageous in terms of computational cost to consider the question of optimal acquisition with respect to the illumination (set of either incident and adjoint
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Illustration of the relationship between the illumination $\phi(x)$ and take-off angle $\varphi(x)$ at target $\mathbf{x}$ in the subsurface. Here we take position $\mathbf{x} = \mathbf{s}$ of a source situated at the surface. Note that the two angles follow different conventions and subsequently they are equal.</figcaption></figure>

directions $\mathbf{p}_s$ and $\mathbf{p}_r$ or illumination angles $\phi_s$ and $\phi_r$). This means that instead of finding directly the optimal positions we find first the optimal illumination and look in a later stage for the geometry on the surface that corresponds to it. This is why an interesting question to ask is the reverse of what has been done so far: for a target in the subsurface, given a certain illumination (whether optimal or not), what is the acquisition on the surface that corresponds to it?

In a homogeneous case and knowing the depth of the target with regard to the acquisition line, the question transforms into a trivial trigonometry problem. However in a heterogeneous case the problem is more complex. One way to tackle this problem is to shoot rays from the target at opposite directions from the optimal illumination and record the positions on the acquisition line at which the rays arrive. In our case we have chosen another approach where we treat the question as a take-off angle problem.

For any source, the take-off angle $\varphi(x)$ is defined to be the shooting angle at the source, which yields the ray linking it to the point $\mathbf{x}$ under consideration (Le Bouteiller 2018). In a 2-D isotropic medium the points along a given ray share the same take-off angle and thus the take-off angles are perpendicular to the direction of the rays. Consequently the gradient of the take-off angle field $\nabla\varphi(x)$ is orthogonal to the gradient of the arrival time $\nabla T(x)$ and the problem can be stated as a simple partial differential equation (Le Bouteiller 2018)

$$ \nabla\varphi(\mathbf{x}).\nabla T(\mathbf{x}) = 0. \quad (10) $$

By considering the target $\mathbf{x}$ as the shooting source and numerically solving eq. (10) we get in the model space the map $\mathbf{x} \mapsto \varphi(\mathbf{x})$ of the take-off angles. Fig. 10 illustrates the relation between take-off $\varphi(\mathbf{x})$ and illumination angles $\phi(\mathbf{x})$. We have chosen the convention for the illumination that the angle is made by the direction of the ray with the downward pointing normal $n_1 = (0, 0, 1)$ whereas the convention used for the take-off is that the angle $\varphi(\mathbf{x})$ is the one the ray makes with the upward pointing normal $n_1 = (0, 0, -1)$; subsequently we have the relationship $\varphi(\mathbf{x}) = \phi(\mathbf{x})$ and solving eq. (10) is equivalent to finding the illumination map $\mathbf{x} \mapsto \varphi(\mathbf{x})$. If we restrict the illumination to the acquisition line (with a fixed depth $z$) we only need to find the horizontal $x$-component of the positions of sources/receivers that correspond to the optimal illumination.

Now that we have discussed the mapping between the acquisition and the wavenumbers in both the homogeneous and heterogeneous cases, it is time to introduce the tools we use in our approach. A corner stone of our method is a concept from computational

geometry called Voronoi tessellation which expresses the problem of regularity of the sampling as a problem of regularity of partition of the sampled domain. Next we review this tool and show how it is relevant to our problem.

## 2.3 Centroidal voronoi tessellation

We want to define a criterion to quantify the repartition spread of the wavenumber point cloud in the wavenumber space. We are going to use Voronoi tessellations to define such metric. While in this study, for the sake of simplicity, we restrict ourselves to the 2-D case, the concepts introduced can be generalized to Euclidean spaces of any dimensions. We explicitly say when the results are restricted to the 2-D or 3-D case.

Let $\Omega \subset \mathbb{R}^d$ be a closed and bounded domain endowed with an Euclidean norm $\|\cdot\|$ and we suppose a constant density function on $\Omega$. A tessellation or tiling is the covering of a surface using geometric shapes. Mathematically speaking (Du et al. 1999) a tessellation of $\Omega$ is the set $\{V_i\}_{i=1}^n$ of subsets of $\Omega$ such that they are disjoint ($V_i \cap V_j = \emptyset, \forall i, j = 1, \dots, n$ such that $i \neq j$) and the domain is the union of their closures ($\cup_{i=1}^n \tilde{V}_i = \Omega$).

A Voronoi tessellation is a type of tiling which is the partition of the domain into regions based on the proximity to a set of points. Let $X = \{x_i\}_{i=1}^n \subset \Omega$ be a set of $n$ points that we call seeds. We define the Voronoi tessellation $\{\Omega_i\}_{i=1}^n$ generated by the seeds $X$ as the partition of $\Omega$ into $n$ Voronoi cells $\Omega_i$ such that each cell is generated by a seed $x_i$ and is defined as the subset of $\Omega$ containing the points closer to $x_i$ than any other seed:

$$ \Omega_i = \{y \in \Omega \mid \|y - x_i\| \le \|y - x_j\|, \forall j \ne i\}. \quad (11) $$

In sampling theory the Voronoi cells correspond to the amount of information each seed represents.

In Fig. 11 we show examples of Voronoi tessellations of a square generated by three different sets of seeds (points). For each Voronoi cell $\Omega_i$, we are interested in the area $m_i$ and the centre of mass $c_i$ (crosses) defined mathematically using surface integrals

$$ m_i = \int_{\Omega_i} dy \quad \text{and} \quad c_i = \frac{\int_{\Omega_i} y dy}{m_i}. \quad (12) $$

In the panel (a) we consider randomly distributed seeds. We note that for an arbitrary choice of seeds it yields arbitrarily sized Voronoi cells. In contrast, in the panel (b) we have examples of regularly distributed seeds and we note that the cells are more uniform in shape and size. The intuition suggests that there exists a link between the regularity of the distribution of seeds and the shapes of cells. The essence of the method we are going to use consists in the iterative update of the positions of the seeds so as to generate a final tessellation where the cells are compact. A byproduct of such regular tessellations is a regular sampling of the domain $\Omega$ by the seeds.

In the literature there exists a misfit function $F: \Omega^n \mapsto \mathbb{R}_+$ which, given seeds $X$, measures the regularity of the generated Voronoi tessellation in terms of cells. It is defined as (Du et al. 1999)

$$ F(X) = \sum_{i=1}^{n} \int_{\Omega_i} \|y - x_i\|^2 dy, \quad (13) $$

and it sums, for each cell, the integral of the squared distance of any point in the cell to the seed. This means that if there exists cells that are larger than the others then the difference in compactness will be amplified by the square distance and thus the
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Examples of Voronoi tessellations. We consider the domain $\Omega$ to be a square in $R^2$ and we put 16 seeds inside. In (a) the seeds are randomly positioned and as a result the tessellation has cells with irregular shapes and sizes. The seeds in (b) are a local minimizer of the energy function $F$ and therefore generate a stable CVT of $\Omega$ where the repartition of the points is regular in spacing. In (c) the tessellation represents an unstable CVT where the seeds correspond to the centroids but they do not minimize the function $F$; as a result the repartition is non-optimal. The Voronoi tessellations have been generated using the SCIPY library (Virtanen et al. 2020).</figcaption></figure>

value of the function will significantly increase. Minimizing $F$ ensures an equilibrium situation where the cells are regular and compact, consequently they are of similar size and each seed $x_i$ is representative of approximately the same amount of information in $\Omega$. An example of uniform tessellation can be seen in Fig. 11(b) where we note that the seeds are regularly distributed and the Voronoi cells are compact. Comparatively to the arbitrary tessellation in Fig. 11 (a), the value of the function is smaller.

The tessellation (b) of Fig. 11 is generated by seeds that are local minimizers of $F$ which explains the regularity of both the tessellation (in terms of sizes, compactness and shapes of the cells) and distribution of the seeds (in terms of spacing of the points). Coincidentally, the seeds and centroids of the associated Voronoi cells share the same positions. Geometrically speaking this corresponds to a special instance of Voronoi tessellation called centroidal voronoi tessellation (CVT): it amounts to finding a placement of the seeds such that they are exactly the centroids of their respective Voronoi tessellations (Choksi & Lu 2020). This is a quite particular situation in the sense that for an arbitrary Voronoi tessellation, the seeds will not coincide with the centroids. To understand the link between minimizers of $F$ and CVTs we need to take a close look at the gradient of $F$. This gradient is defined in terms of volumes and centroids of the Voronoi cells through the expression

$$ \nabla F(X) = [2m_i (x_i - c_i)]_{i=1}^n. \quad (14) $$

The gradient vanishes if $x_i = c_i$, for all $i = 1, \dots, n$. So a CVT can be analytically described as the tessellation generated by seeds $X^*$ that are a critical point of the misfit function $F$. Note that a critical point does not guarantee the compactness of Voronoi cells because it can simply be a saddle point. This is why a minimizer of $F$ necessarily generates a CVT but the reverse is not true: a CVT is not necessarily associated with a minimizer of $F$. For this reason Liu et al. (2009) distinguishes between stable CVTs which are local minima of $F$ and unstable CVTs.

An example of unstable CVT is shown in Fig. 11(c) where we see that the seeds do coincide with the centroids but the rectangular

cells are not compact. As a consequence, the gradient of the energy function is indeed null but the value of the function is not minimal. The effect is an irregular sampling of the domain where the sample points are not optimally distributed: the horizontal sampling is sparser than the vertical sampling.

Moreover optimal CVTs are defined as Voronoi tessellations generated by the global minimizer of $F$. However the CVT problem is not convex and for our application any local minimizer will be enough as it offers a steady state tessellation. For example the tessellation in Fig. 11(b) is not the optimal CVT but is still adequate in terms of regularity of the sampling. Consequently, we are interested in stable CVTs to get a regular sampling of the wavenumber domain $\mathcal{W}$.

To minimize $F$, we will use a local optimization method. Due to the piecewise and complex nature of the function $F$, it was for a long time believed that it is non-smooth and this discouraged the use of Newton-like optimization methods. However Liu et al. (2009) proved in 2-D and 3-D that if the domain $\Omega$ is convex then the function $F$ is of class $C^2$. This new result has prompted the use of Newton methods for finding stable CVTs. However, Newton method is not recommended for large scale problems because of the computational cost of the Hessian. If we link the CVT to our acquisition problem then we need to optimize the sampling for $n = N_s \times N_r$ points which can increase rapidly by just adding a source or receiver. This is why we rely on quasi-Newton methods, more precisely the I-BFGS method, which is the limited-memory version of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm (Nocedal & Wright 2006).

We take advantage of the CVT characterization of the sampling problem to improve the regularity of the wavenumber cloud. In Fig. 12 we consider different samplings of the wavenumber domain and their effect on the generated Voronoi tessellations. The wavenumber domain $\mathcal{W}$ is the same as in Fig. 7 in which we compared the regular and optimized acquisitions. On the top panel, we show the wavenumber cloud associated to the regular acquisition. Again the sampling is irregularly distributed inside the envelope and the effect on the tessellation is Voronoi cells of varying sizes.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Voronoi tessellations of the wavenumber domain $W$. We reuse the same dimensions as in Fig. 7 (number of sampling points, envelope) and compare two samplings in terms of tessellation: one generated by the regular acquisition of Fig. 7 (top) and the other one a local minimizer of $F$ (bottom). The distribution of the points inside the envelope is more uniform in the second case and on the panel (b) where we plot the Voronoi tessellations of $W$ this corresponds to a CVT. So by minimizing the value of $F$ we manage to attain a regular sampling of the wavenumber space.</figcaption></figure>

On the bottom panel, we have in contrast a stable CVT of $W$ generated by a local minimizer of $F$. We achieve to halve the value of $F$ and consequently to gain both regularity of the sampling and tessellation. We remind the reader that here we optimized the sampling without the source-receiver formulation (2) so we considered $N_s \times N_r$ degrees of freedom. Therefore the sampling is ideal and not representative of any fixed-spread acquisition. This is another reason why we do not expect to attain optimal CVTs: the global minimum of $F$ is unreachable for our problem because of the limitations imposed by the interdependence of the wavenumbers in the source-receiver formulation.

In light of this development, we present in the remainder of this section our approach for the optimal design of acquisition.

## 2.4 Approach

The approach that we propose takes advantage of the analytical characterization of the CVT problem. The idea is to combine the optimal acquisition problem with the CVT problem to devise a new objective function which provides a measure of quality to an acquisition with regard to the regularity of the wavenumber sampling.

As a prior knowledge we set the position of the target, the characteristics of the medium (dimensions, velocity), the reference frequency, the limits of the deployment and the number of sources and receivers. Again we deal with fixed-spread acquisitions, so for $N_s$ sources and $N_r$ receivers we get $N_s \times N_r$ wavenumbers. However,

because we optimize the positions of the sources/receivers, we only have $N_s + N_r$ degrees of freedom.

Given an acquisition $(S, R)$, we compute the associated wavenumber content $K$ through the map $G$ defined in (3). We deal with surface acquisitions so the depth $z$ will be set to zero. We can choose to represent the acquisition either by the positions $x_s$ and $x_r$ of the sources and receivers or by its illumination angles. For the 2-D case we have decided to optimize the acquisition in terms of illumination angles $S = (\phi_s)_{s=1,...,N_s}$ and $R = (\phi_r)_{r=1,...,N_r}$. As explained later, this makes it possible to confine the use of ray-tracing and take-off angles solvers to a pre- and post-processing steps and thus reduce the computational cost of the procedure.

From here we define a new function $F \circ G$ which will be the composition of the CVT energy function $F$ and wavenumber map $G$:

$$ F \circ G: \begin{array}{c} D_s^{N_s} \times D_r^{N_r} \\ (S, R) \end{array} \rightarrow \begin{array}{c} W^{N_s \times N_r} \\ K \end{array} \rightarrow \begin{array}{c} R_+ \\ F(K) \end{array}. \quad (15) $$

Here $D_s = D_r = [\phi_{min}, \phi_{max}]$ are the domains of all possible angle values and denote in fact the illumination range. For an acquisition $(S, R)$, the function assesses the quality of the geometry with regard to the regularity of the wavenumber cloud $K$. Indeed, it computes $F(K)$ and therefore links the sampling to the regularity of the corresponding Voronoi tessellation. Thus minimizing $F \circ G$ is equivalent to minimizing the CVT of the wavenumbers in the domain $W$. It means that we have a function that takes as input the acquisition
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Illustration of our approach. Knowing the position of the target we compute the wavenumber content associated with the acquisition using the map G. Knowing also the maximum offset we are able to draw the envelope of the wavenumber cloud and from here use the points inside as input to the CVT energy function F which assesses the regularity of the Voronoi tessellation, thus also the regularity of the sampling inside the envelope. This workflow gives a quality measure for the acquisition geometry with regard to our regularity criterion.</figcaption></figure>

and outputs its fitness with regard to our chosen criterion. Fig. 13 illustrates this workflow.

Finally we express the problem of finding an optimal acquisition as the minimization problem

$$ \min_{(S,R)} F \circ G(S, R). \qquad (16) $$

We have decided to use a local optimization strategy to solve (16). Because we have settled for an angular characterization of the wavenumbers, the map G is expressed through eq. (4) using trigonometric functions and is therefore of class $C^\infty$. However the wavenumber domain W defined in eq. (6) is not convex as it can be seen in the shape of the envelope in Fig. 6. So F is not guaranteed to be $C^2$. Nonetheless Liu et al. (2009) argue that in non-convex domains it is still $C^2$ in most situations encountered in optimization. Indeed, F is not smooth in the situation where a face of Voronoi cell is parallel to a segment of the envelope. But with a large enough number of seeds this situation is unlikely to happen. We can therefore assume that $FoG$ is $C^2$ almost everywhere and similarly to the CVT problem we choose to use the same quasi-Newton strategy as for the CVT, namely the I-BFGS method.

The gradient of F is known through formula (13) and because G is expressed with simple trigonometric functions then its Jacobian is easily derived. So we can express the gradient of $F \circ G$ using the chain rule

$$ \nabla(F \circ G) = J_G \cdot \nabla F(G). \qquad (17) $$

The main incentive for the choice to solve (16) with regard to the illumination instead of positions lies in the reduction of the computational cost. In fact, when using the positions of sources and receivers, extra calculations are needed at each iteration of the optimization process to get the $p_s$ and $p_r$ vectors in eq. (2). In a homogeneous medium those calculations are trivial. However for an arbitrary heterogeneous medium, as seen in subsection 2.2, we need to numerically solve an eikonal equation and bear the underlying computational cost at each iteration since the positions of sources and receivers would have been updated. An alternative would be to compute the take-off angle map once by numerically solving eq. (10), and reduce the finding of the $p_s$ and $p_r$ vectors as a simple computation using the illumination angles that we can get easily from the take-off angles map. For practical reasons we have opted for the use of ray-tracing to estimate the illumination directions. But with our angular choice, ray-tracing needs to be performed only as a pre-processing step to compute the illumination from an initial acquisition. Even this step can be skipped if we decide to use an initialization directly defined with angles. Similarly the take-off

angle problem needs to be solved as a post-processing to get from the optimal angles $\phi_s$ and $\phi_r$ the corresponding positions of sources and receivers.

## 2.5 Implementation

Given the relatively low computational requirements of the 2-D case, compared to 3-D, we have decided to implement the method in Python using NUMPY (Harris et al. 2020) and other tools. The wavenumber map G takes as input an array of angular values representing illumination and is expressed using trigonometric functions so its implementation is straightforward. To compute the CVT energy function F and its gradient, knowing the Voronoi tessellation of the domain, we need to evaluate surface integrals that appear in eqs (13) and (12). To this end we use the divergence theorem to transform the surface integrals over polygonal cells into simpler line integrals along the edges.

Concerning Voronoi tessellations there are many methods to generate them. There exists a link between them and Delaunay triangulations: they are dual. The vertices of the tessellation are the circumcentres of the triangles in the Delaunay triangulation. Inversely the vertices of the Delaunay triangulations are the input points of the Voronoi tessellation. So if we manage to get the Delaunay triangulation we can get the Voronoi tilling.

To compute the Voronoi tessellations we use SCIPY (Virtanen et al. 2020). It uses the Qhull library in C++ (Barber et al. 1996). This latter implements the Quick Hull algorithm for the computation of convex hulls of a set of points. This algorithm is of worst case time complexity $O(n^2)$ with n being the number of points; in our acquisition problem we would have $n = N_s \times N_r$. Qhull computes the Voronoi tessellation of a set of points from their Delaunay triangulation. To compute the triangulation, Qhull uses the relation that exists between Delaunay triangulations and convex hulls: The points in $R^2$ are lifted in $R^3$ by a height $z = x^2 + y^2$ to a paraboloid, and the program computes the convex hull of the set of points in $R^3$. Afterwards it projects the hull into the plane and this creates the Voronoi tessellation (Chen et al. 2022). SCIPY only provides tessellations in an unbounded space (with edges going to infinity) and we have to clip the vertices to restrict it to our domain. To manipulate such geometric objects in 2-D we have used the Shapely package (Gillies et al. 2024).

For the ray-tracing we have chosen to use the Hamiltonian fast marching library (Mirebeau & Portegies 2018), a fast-marching method to solve eikonal-type equations. The method is implemented in C++ with an interface in Python. It takes as input the
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Two models associated with the targets—they are composed by a homogeneous background with a Gaussian anomaly at the position of the target.</figcaption></figure>

velocity model, characteristics of a box-shaped grid (number of points, grid dimensions, step size) and the position of the source. It then computes the arrival times $T$ and its gradient $\nabla T$ at the grid points. To extend the values outside the grid points we again use interpolation.

After the minimization process (16) we use the FSM-DG program developed by Le Bouteiller (2018) to solve the take-off angle eq. (10) by considering the source to be the target $x$ for which we optimized the acquisition. This computes the values of the take-off angle $\varphi$ at points of a grid covering the whole model.

For performing FWI and computing FWI gradients, we use TOYxDAC_TIME, a finite-difference based full waveform modelling and inversion code developed in the SEISCOPE project (Yang et al. 2018). The solver used for the optimization is the SEISCOPE optimization toolbox (Métivier & Brossier 2016) which implements the I-BFGS method in Fortran.

In the next section we illustrate our work with numerical examples.

## 3 NUMERICAL EXPERIMENTS

In this section, we present numerical experiments which illustrate the scope of application of our method. For each application, we consider one regular and one optimized acquisition in order to compare them in terms of acquisition geometry, wavenumber content (both sampling and Voronoi tessellations) and imaging result (first FWI gradient and inversion).

Throughout those applications, the acquisitions are deployed at the surface and consist of 16 sources and 17 receivers. The subsurface model is 17 × 3.5 km and increasingly varies in complexity. At first we consider a homogeneous model, then move to a linear vertical velocity gradient and finally consider the Marmousi model (Martin et al. 2006), in order to represent a realistic set up. For each model, we consider 2 or 3 target points located at different spatial positions and optimize the acquisition for each of them.

### 3.1 Homogeneous background models

We start with a simple homogeneous background model with a constant velocity 1500 m s⁻¹. Two targets are considered, located at the centre (with coordinates (8.5, 1.5) km) and right bottom corner (with coordinates (13, 3) km), respectively. For each target, the model consists in a homogeneous background with a small Gaussian-shape (diameter ≈ 300 m) velocity anomaly (peak velocity of 4000 m s⁻¹) at their position (Fig. 14).

We optimize the acquisition for each of the target and run a FWI for 50 iterations to retrieve the anomaly from the homogeneous background.

In Fig. 15 we compare a regularly spaced acquisition to the optimized acquisition for the central target. The geometries and wavenumber samplings are similar to the ones shown in Fig. 7 where we compared acquisitions in a homogeneous medium. The optimized geometry, in an attempt to regularize the illumination, puts more sources and receivers density above the target and then sparsifies the density with distance. Because of the difference in offset range (previously the range of deployment was 4 km for Fig. 7 and here it is 17 km), the wavenumber envelope is larger and the irregularity of the initial wavenumber sampling is more pronounced. The irregularity is visible in the initial Voronoi tessellation, where the cells are of various sizes. In contrast, the optimized wavenumber sampling and tessellation are regular. Note how the energy function decreased almost by a factor of five, from a value of 5.7 down to 1.2 after optimization. We then compare the imaging results with the first FWI gradients. We can see that for the regular acquisition, the gradient is characterized by an update which is more spread out than for the optimized acquisition and displays stronger artefacts in the reconstruction, particularly visible on the horizontal profile. Those artefacts in the gradient are transferred to the inversion result after iterations. In the last row of Fig. 15, the optimized acquisition provides a smoother focused update at the target and thus a better reconstruction of the anomaly. Regardless of the optimization, we still reconstruct a Cardinal-Sine-like shape on the vertical section due to the limited illumination and wavenumber coverage.

In Fig. 16 we show the results for the bottom-right target. Almost the same conclusion can be drawn: the optimized geometry has denser sources and receivers above the target and a more regular wavenumber distribution inside the envelope, corresponding to a decrease in the energy function, and a more focused reconstruction of the anomaly.

We can still point out some interesting differences. For a deeper target the initial acquisition produces a less irregular sampling compared to the shallow case. This smaller difference between the initial and optimal values of $F$ translates into a lesser improvement in the imaging. Note also that the optimal geometry is shifted to the right of the deployment range and this to be able to sample the illumination range [$\phi_{min}, \phi_{max}$] = [−54°, 77°] which is not symmetric in relation to the origin; for comparison, in the centred case the illumination range was [−80°, 80°]. This asymmetry has the effect of producing a wavenumber envelope which is tilted to the right and the initial wavenumber cloud concentrates the wavenumbers against the rightmost circular arc of the envelope, given that there are more sources and receivers on the left of the target and therefore
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Results for the homogeneous case with a target in the centre. The left column is dedicated to the regular acquisition and the right one to the optimized acquisition. The geometries are shown on the first row. Below, the corresponding sampling and Voronoi tessellations of the wavenumber domain $\mathcal{W}$. The next rows show the first FWI gradients (L2-normalized) and the perturbation (final model minus starting model) after 50 iterations of FWI. We add zoomed images around the target for clarity.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Results for the homogeneous case and target in the bottom right. Compare with Fig. 15, where the target is centred.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 17.** 1-D velocity models with two targets at different positions—at the positions of the targets are Gaussian anomalies. They will be used as the true model for the inversion while the 1-D background will be the initial one.</figcaption></figure>

a larger sampling of the right half of the wavenumber space. The equal repartition of the sources/receivers on each side of the target in the optimized acquisition indicates the fitness of the geometry with respect to the target and produces wavenumbers that are regularly distributed inside the tilted envelope. Moreover, the tilt of the wavenumber sampling has the effect of producing tilted reconstructions in the gradient and inversion images, no matter if the acquisition is optimized or not. Another contrast with the centred case worth mentioning is the smaller size of the wavenumber envelope. This difference comes from the variation in depth of the target: previously we had a shallower target so the illumination range was larger.

## 3.2 1-D velocity gradient background model

In this second application, we consider a 1-D velocity gradient model with velocity increasing from 1500 to 4000 m s⁻¹. For each target, we define the true model to be the 1-D background on which is added the Gaussian anomaly at the target positions, as shown in Fig. 17. In a similar way as the homogeneous case, we compare the acquisitions and FWI results at iteration 50 in Figs 18 and 19.

The optimized geometry resembles the one optimized for a centred target in a homogeneous medium, shown in Fig. 15, in that the spacing between sources and receivers progressively increases with the distance to the target, although the accumulation is less pronounced than the homogeneous case. Despite the heterogeneity, a certain symmetry is kept due to the invariance of the velocity model in the x-direction.

Similarly to Fig. 8 where we show the effects of the velocity heterogeneity on the wavenumber cloud, the wavenumber envelope is characterized by this particular shape where the side circular arcs of the envelope cross each other. Compared with the homogeneous case, the covered area is larger due to the wider range of illumination. In a vertically increasing velocity model such as this one and a wide range of offset on the surface, we can expect significant undershooting of the target. This is confirmed by the significant proportion of negative $k_z$ wavenumbers in the top-half of the wavenumber space. Those negative wavenumbers are associated to a large contribution of illumination directions $\mathbf{p}_s$ and $\mathbf{p}_r$ pointing upwards. You can note that this avoids the cardinal-sine effect (oscillations along $z$). For the optimized acquisition the sampling is, as expected, more uniform and we achieve a threefold decrease of the CVT energy function. When comparing the FWI gradients, we achieve a stronger update with the optimized acquisition and in turn the inversion image shows a larger and more focused anomaly.

In Fig. 19 we show the results for the bottom right target. Similarly to the homogeneous case, the target is deeper and uncentred, creating

a smaller and tilted envelope. Because the target is close to the limits of the model, we pay attention to consider dimensions of the model such that the whole possible paths are not hampered by the model limits. Therefore, the model we considered is larger than the one displayed in the figures. For this case again, the use of the optimized acquisition allows to halve the value of the objective function and get a reconstruction of the anomaly which is more focused and with less artefacts.

## 3.3 Marmousi model

In this third test, the Marmousi II velocity model is considered. Marmousi is a pioneering synthetic model that has been realistically enriched in its second version by Martin *et al.* (2006) to contain several hydrocarbon reservoirs trapped in sand layers at different depths with varying sizes and structural complexities. In order to test our approach in this realistic context, three of those hydrocarbon reservoirs are selected as targets. The positions of the three targets are located at positions $T_1 = (10.8, 2.85)$ km, $T_2 = (3, 1.1)$ km and $T_3 = (10.45, 1.13)$ km in the model. Fig. 20 shows the true Marmousi model, a smoothed version and the three targets. The smoothed model represents a low-resolution prior knowledge of the medium, used for the optimization of the acquisition and as starting model for FWI.

In Fig. 21 we display the results of the inversions at iteration 20. As for the two previous test-cases, the optimized acquisition shows a gathering of sources and receivers above the targets. It is interesting to note in this realistic set-up, that the general model reconstruction is better for the regular acquisitions. However, the three optimized acquisitions provide improved local reconstruction around the targets. For the target $T_1$, the imaging using optimized acquisition provides stronger details to the layering. For the targets $T_2$ and $T_3$, the reconstruction appears less noisy and more focused with the optimized acquisition.

To highlight those differences, we show in Fig. 22 the error between the inversion and the true model. We also compute the $L^2$ norm of this error in the zoomed-in region to assess the local improvements. We can see for each target that we locally reduce the error around the target when we use optimized acquisitions, either by gaining stronger updates (targets 1 and 3) or by reducing the artefacts (targets 2 and 3). This reduction in error value confirms the better local imaging of the medium with the optimized acquisitions.

# 4 DISCUSSION

After presenting the theory behind our method and three test-cases of increasing complexity, we would like to discuss several interesting points related to our approach.
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Results for the 1-D case with a target in the centre. The results of the regular acquisition are compiled in the left column while the results of the optimized acquisition are in the right column. The geometries of acquisition are shown in the first row while the sampling and tessellation of the wavenumber space are shown in the second row. On the third and fifth row we show, respectively, the first FWI gradient (L2-normalized) and perturbation (final model minus starting model) at iteration 50. Between them we show, for the sake of clarity, zoomed in images around the target, corresponding to the square regions.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 19. Results for the 1-D case—target in the bottom left corner. Compare with the centred case in Fig. 18.</figcaption></figure>

## 4.1 Scope of the method

The numerical results presented in the previous section shows the successful implementation and use of our method to find optimal acquisition geometries for the improvement of the imaging of a

target point in the subsurface. Throughout the section we have so far shown results with few sources and receivers, namely 16 sources and 17 receivers. In all examples, the optimized acquisitions provided significant improvement over the regular geometry. We conduct
<figure><img src="image_21.png" /><figcaption>**Figure 20.** The true Marmousi II model (top) and the initial smoothed model (below) used for the FWI, along with the three targets located at different hydrocarbon layers.</figcaption></figure>

here another experiment, using the same Marmousi settings than in Section 3.3, except we consider now more sources and receivers in the acquisition. Instead of 16 sources and 17 receivers, we now consider 128 sources and 129 receivers.

As previously, the results are compiled in Figs 23 and 24. We succeed in optimizing the acquisition. Even if the large number of sources (downward-pointing triangles) and receivers (upward-pointing triangles) makes them difficult to be distinguished, we can still differentiate between the regular and the optimized acquisitions since optimized geometries display a distinct concentration above the target while the regular geometry occupies the whole surface.

For imaging, the dense regular acquisition already provides a very good global resolution of the model. We are able to improve the reconstructions for targets 1 and 3 with an optimized acquisition, as can be seen in the reduction of the local error value. However, for the second target, the optimized acquisition fails at providing a gain compared to the regular case (error of 1.2 compared to 0.71). More importantly, the zoomed-in image around target 2 shows stronger artefacts with the optimized acquisition. We suspect that this results from the combination of two facts :

(1) the shallow depth of this target makes possible significant undershooting from reflections coming from the strong velocity layer underneath (Mora 1989),

(2) the optimized acquisitions which do not take benefit of the whole offset range on the surface.

Thanks to this undershooting, the range of sampled wavenumbers is much richer than the one we can expect from the surface only. However, the high frequency approximation of the wave equation we consider in the eikonal, used for estimating the range of sampled wavenumber, does not take into account the undershooting associated to second-order scattering effects. Consequently, our method operates under pessimistic assumptions to define the wavenumber envelope and the effect of the reflections is not taken into account in the optimization of the acquisition.

If we consider the improvement of the imaging in the two other targets, it has to be noted that the gain from the optimization is clearly weaker compared to the sparser acquisition case. This means that our method does not provide enough gain for dense acquisitions and its advantages lie in scenarios where the number of sources and receivers is limited and the global goal is to maximize the effect of each device to image the target.

Given the presented workflow of our target-oriented approach, in which a prior knowledge of both the position of the target of interest and the velocity model (to estimate the wavenumber sampling) is

necessary, our method appears to be well-suited for monitoring applications.

## 4.2 Effect of the initialization

Given the nature of our local optimization method for the acquisition, a natural question is the effect of the initialization on the convergence towards an optimal solution. To tackle this question we consider an experiment where we optimize the acquisition for one target and by using different initializations. Similarly to Fig. 7, we consider an acquisition of 10 sources and 11 receivers on a deployment range of 4 km at the surface. The model is homogeneous and the target is at the centre of the range and at depth 1 km. Since we can express the wavenumbers in terms of either positions or angles, we can choose the initialization also by expressing it using positions or angles. We chose two regular acquisitions, one in terms of positions and the other in terms of illumination. We compare them with random initializations. We generate the random acquisitions in two ways: first the positions follow a uniform distribution on the range [0, 4] of the possible values of the x-component; for the second, the angles are generated uniformly on the range [−63°, 63°] representing the illumination range.

For each type of random initialization, we generate 49 realizations of starting points. In Fig. 25, we show a plot of the values of the objective function *F*: at the initialization (a) and after optimization (b). For comparison, we plot as well the values of the function for regular starting points (crosses). Each realization is plotted as a dot (blue for random positions and green for angles) and we also show a barplot of the data (each of the four sections represent a quartile of the realizations). When looking at the initial values in Fig. 25(b), we note that all the random realizations start with a worse energy function than the regular initializations. After optimization, this tendency is reversed for both the position and angle cases since 75 per cent the random realizations yield a value of *F* lower than regular initializations, as it can be seen in Fig. 25(b), where the crosses are positioned just below the third quartile. This hints at the role the initialization plays in the method: an initialization which already provides a good sampling would have the effect of generating an optimized acquisition which gravitates around a local minimum not far away from the initialization.

For illustration, we compile some acquisitions of this experiment in Fig. 26. We show four initial acquisitions: the regularly spaced, the one with the regular illumination, and the two geometries which yield the median *F* value after optimization, among, respectively, the random positions and random angles. For each case we plot both the initial and optimal geometries obtained (top two panels), along with their *F* values. In the two bottom panels we show the corresponding wavenumber samplings.
Results inversion for target 1 (x=10.80 km, z=2.85 km)

<figure><img src="image_22.png" /></figure>

Results inversion for target 2 (x=3.00 km, z=1.10 km)

<figure><img src="image_23.png" /></figure>

Results inversion for target 3 (x=10.45 km, z=1.13 km)

<figure><img src="image_24.png" /><figcaption>**Figure 21.** Inversion results for Marmousi with 16 sources. For each of the three targets we compare the inversion results after 20 iterations (taking as initial model the smoothed version of Fig. 20) of a regular and an optimized acquisition. We zoom inside the rectangles to highlight the local improvements around the targets.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure 22.** Errors of the inversion for Marmousi with 16 sources. For each target we compare the errors (true model minus final model) of the inversions after 20 iterations. We zoom inside the rectangle region around the target and also compute the L2-norm of the errors in this region in order to highlight and assess local differences.</figcaption></figure>
Results inversion for target 1 (x=10.80 km, z=2.85 km)

<figure><img src="image_26.png" /></figure>

Results inversion for target 2 (x=3.00 km, z=1.10 km)

<figure><img src="image_27.png" /></figure>

Results inversion for target 3 (x=10.45 km, z=1.13 km)

<figure><img src="image_28.png" /><figcaption>**Figure 23.** Inversion results for Marmousi with 128 sources. Compare with Fig. 21 where we had 16 sources.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>**Figure 24.** Errors of the inversion for Marmousi with 128 sources. Compare with Fig. 22 where we had 16 sources.</figcaption></figure>

On the top row we have an initial acquisition that is regular with
respect to positions. The sources and receivers are alternating and
regularly spaced. For the regular-illumination acquisition, we can
see that the spacing is no more regular but rather increasing the

further we go from the target. We note that the initial *F* value is
smaller for the regular illumination than the regularly spaced one.
However, after optimization, both yield a similar energy function,
so the optimized two acquisition geometries are almost identical.
<figure><img src="image_30.png" /><figcaption>**Figure 25.** Two barplots of the *F* values of random initializations: (a) at the initialization and (b) after the optimization. We plot the values of *F* for each random realization (dots) and compare it with the regular initializations (crosses). Each section of the barplot represents a quartile of the data.</figcaption></figure>

For the random initializations, we can see that we have arbitrary geometries in both the position and angle cases. This has the effect of generating wavenumber samplings of poor quality with many holes inside the envelope and comparatively high values of the initial *F* function. Despite this unfavourable starting point, after optimization we are able to recover geometries that are optimal in terms of sampling, such that we get objective functions which are even weaker than the two previous optimal acquisitions. The sampling is good at the centre of the envelope (less spacing between the points), although it still suffers from accumulation at the bottom arc.

We learn from this study that despite the non-convexity of the problem, the method converges most of the time towards satisfying results, independently of the quality of the initial start. This can be seen statistically in the important decrease of the objective function in all the realizations in Fig. 25 and graphically from the example acquisitions of Fig. 26, where the initial wavenumber samplings were far from satisfactory. Moreover, we note that the largest reductions in the objective function were achieved using random initializations. To gain the lowest value after optimization, a strategy is to run the method for different initializations and take the result which provides the best outcome. This strategy is reasonable when working on small enough acquisitions. When dealing with large acquisitions, a pragmatic view that we hold is to use regular initializations in terms of angles in order to ensure good quality optimized acquisitions, because this would represent already a good starting point. In our workflow, we have also decided to express the acquisition in terms of angles, and choosing an initialization which is uniformly distributed in illumination will have the advantage of skipping the pre-processing step of linking an initial acquisition (expressed in positions) to the illumination, thus avoiding the ray tracing calculations. We do not recommend using initial acquisitions that are regular in positions, because the initial wavenumber cloud is not always good in terms of distribution, especially for large offsets, and the heterogeneity might affect strongly the initial sampling in an unpredictable fashion.

## 4.3 Stability with respect to target position

Another question worth asking concerns the stability of the method with respect to parameters of the problem. Specifically we discuss here the stability of the results with regards to the target: what is the effect of its position on the workflow and the optimal solution?

Does a small change in position results in a drastic modification of the optimal acquisition geometry?

We conduct a numerical experiment and use the same configuration as before: an acquisition of 10 sources and 11 receivers on the surface of a 4 × 2 km homogeneous model. We consider 10 different positions of the target. We have a diagonally moving position from the top left corner of the model to the bottom right, and we optimize the acquisition for each of them. For the first top left target we use a regularly-spaced acquisition as initialization and the next targets will be slightly incremented to move diagonally. To preserve a continuity between the geometries we use as initialization the optimized geometry for the previous target.

In Fig. 27 we show the optimized acquisitions for each target, while in Fig. 28 we show the optimized wavenumber contents and the envelope associated to the targets. The first target is shallow and uncentred, so the optimized acquisition is centred on the left side to be aligned above the target. In the wavenumber cloud we can note that the envelope is tilted and despite this, thanks to the optimization, we are still able to regularly sample the interior of the envelope.

We can see that for each spatial increment of target position, there is a small incremental change of orientation to the wavenumber cloud. This is due to the change in tilt of the illumination. Moreover, the size of the envelope is decreasing with the increase of the distance of the target from the surface. This is because the offset in the acquisition remains constant while the depth of the target increases, and this affects the illumination range, which becomes narrower with depth. Consequently, the size of the envelope decreases, as has been discussed in the geometric characterization of the wavenumber cloud.

Importantly, we note that an incremental change to the target produces a gradual update to the optimized acquisition. Therefore, a small change in the parameters of the problem leads to a small change of the solution. This is also a sign of the local nature of the optimization, and again taking an initialization which is relatively good leads to few changes to the acquisition and a convergence towards a local minimum.

## 4.4 Mathematical perspective

We discuss here the method from its mathematical and implementation perspectives.
<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /><figcaption>**Figure 26.** Plots of some acquisitions of the initialization experiment. For each initial geometry, we present both the optimized acquisition after the optimization process and the corresponding wavenumber sampling.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>**Figure 27.** Acquisition geometries optimized for different targets. The characteristics of the medium and acquisition are the same, and the difference is the diagonally moving target. Each acquisition is optimized using as initialization the optimized geometry of the previous target (from the top left to the bottom right).</figcaption></figure>

The objective function $F \circ G$ is analytically defined; the CVT energy function $F$ can be expressed as line integrals along the edges of the Voronoi cells. As discussed earlier, we have chosen to optimize for the illumination angles instead of positions, so we express $G$ with simple trigonometric functions. Additionally, expressing the

map via the positions formulation would also use analytical functions when the medium is either homogeneous or the velocity model is known and expressed by equations for which we can find an analytical solution for the ray tracing. We need to deal with numerical solutions when considering heterogeneity and the stability of the
<figure><img src="image_36.png" /><figcaption>**Figure 28.** Wavenumber samplings of the acquisitions in Fig. 27.</figcaption></figure>

optimized acquisitions depends on the quality of the tools used in the pre- and post-processing phases of the workflow.

Moreover, we concluded that the objective function is $C^2$ almost everywhere. Because the wavenumbers are dependent on the geometry of acquisition, we get duplicated seeds if two couples source/receiver generate the same wavenumber. This is due to the

symmetry of the function $k(s, r) = k(r, s)$ in eq. (2), so the singularity happens when there exists two couples where the positions of the sources and receivers are swapped. To remedy this singularity problem, we decide to test if wavenumbers values are duplicated by simply comparing the points two by two. In case of duplicated wavenumber point, a practical solution is to randomly
<figure><img src="image_37.png" /><figcaption>**Figure 29.** Optimization of the wavenumber sampling of the acquisition when considering a larger envelope (here a disc). The two plots show the wavenumber sampling of the initial acquisition and the 'optimized' one. We can see that the method failed here since it did not change much of the initial sampling so it did not optimize the distribution. Compare this to the result of Fig. 7.</figcaption></figure>

perturb the duplicates by a small amount. Because we have so far worked with small acquisitions, the algorithmic complexity of this search has not been a concern. If we would consider denser acquisitions, the computational cost of these operations could escalate rapidly. Therefore, a better strategy would need to be applied to either manage the singularity or simply detect and remove the duplicates.

Another undefined behaviour occurs when the wavenumbers are outside the domain $\mathcal{W}$. This can happen when the envelope is not representative of the boundaries of the domain (for example if the illumination range [$\phi_{min}, \phi_{max}$] has not been estimated correctly) or when it is too tight (so it does not account for numerical variations in the seeds or the small random perturbations that might apply on them).

To circumvent this difficulty, we might be tempted to define an envelope which is larger than the circular arcs one. This solves the problem of out of bounds wavenumbers but in our experience this creates another problem with the convergence of the optimization. We illustrate this with a small experiment where we want to optimize the regular acquisition of Fig. 7. In this example, we optimize the regularity of the wavenumbers inside of the whole disc $\bar{D}(O, 2k_0)$ of centre the origin and radius $2k_0$ instead of the tighter domain $\mathcal{W}$. In Fig. 29, we show the Voronoi tessellations of the disc with both the initial wavenumber sampling and the optimized one. We can see that practically no change has been applied to the sampling except for the two side arcs which means that the acquisition has not been truly optimized. Here, the 1-BFGS optimization algorithm actually stops at the second iteration after a line search failure. We interpret this behaviour by the fact that in a loose envelope the seeds at the boundary of the actual domain $\mathcal{W}$ have the most impact on the CVT energy function, and the optimization prioritize them over the other seeds and most probably terminate at a local minimum. This highlights the importance of an adequate definition of the wavenumber domain $\mathcal{W}$. In our implementation, a marginally larger envelope has been used to alleviate possible problems of the wavenumbers going out of bounds.

## 4.5 Possible developments and path towards 3-D

Many interesting developments in the optimization of the acquisition have been left out in this work and are the subject of a companion paper extending our work.

For example, we have not taken into account deployment constraints on the geometry of the sources and receivers. Indeed, when

designing seismic surveys we are faced with restrictions such as inaccessible zones, either because of physical obstacles (infrastructure, geographic features) or legal ones (protected natural habitats). A way to add those limitations into the method is to incorporate them as constraints and penalty functions in the minimization problem (16). However, it is more intuitive to define those constraints on the positions of the sources and receivers rather than the illumination, therefore a redefinition of the map $G$ in terms of positions is required.

Despite the target-oriented nature of the method, it is still possible to consider a multitarget approach. For this we optimize the acquisition for several target points by defining the objective function as a sum of functions $F \circ G$, each one associated to a target. This is trivial in terms of implementation and can be easily parallelized due to the independence of the functions. This opens up the possibility to consider targets as volumes instead of just discrete points: to evaluate the fitness of the acquisition over the volume we need to take the integral of the objective function over said volume. Using numerical integration schemes, we can discretize the volume integral into a weighted sum and apply the multitarget approach.

However, the current strategy of optimizing the acquisition based on illumination angles at individual targets does not allow for geometry optimization across multiple targets. This is because the input values in the method are associated to specific subsurface positions. Again, to extend the approach for multiple targets, it is necessary to adopt a new strategy that optimizes the acquisition based on target positions instead. In the 3-D extension, we go beyond the illumination angles parametrization and propose a new strategy to optimize the acquisition in terms of positions, and this in order to avoid the limitations imposed by this choice of strategy (deployment constraints, multitarget).

Most importantly, in this work we have exclusively worked on the 2-D case in order to develop and illustrate our method more easily. There is however no barriers to the extension to the 3-D case, which has become over the years the standard in seismic acquisition. This would lead to interesting developments and geometries when working on a 2-D plane, as opposed to working on a simple line.

All of those developments and more are presented in a companion paper dedicated to the extension of our method to the 3-D case. In the 3-D extension, we will go beyond the illumination angles and propose a new strategy to optimize the acquisition in terms of positions, and this in order to avoid the limitations imposed by this choice of strategy (deployment constraints, multitarget).
## 5 CONCLUSION

We propose a new method for the optimal design of acquisition geometry, with an imaging criterion based on the improvement of the wavenumber sampling at a target point. This sampling is estimated using the theory of diffraction tomography. The method describes the sampling problem as an instance of Centroidal Voronoi Tesselation. The workflow describes a new objective function which, when minimized, provides geometries that are optimal with regard to our regular criterion.

From the numerical perspectives, we have succeeded in implementing our method in the 2-D case: given characteristics of the problem (position of the target, deployment range, prior-knowledge about the velocity model), we are able to find acquisition geometries that are optimized for the target. We have numerically validated our method in the context of FWI with different synthetic tests, where we compare optimized and regular acquisitions in terms of wavenumber sampling, value of the objective function and local quality of the imaging.

In all our tests we achieve to regularize the sampling, and as a result to improve the imaging at the target. However, when considering dense acquisitions the regular geometries provide an already good global resolution so the gain from the optimization is weaker. Moreover, for one test we had a worse resolution of the target. We attribute this to the fact that our method only considers, in a pessimistic manner, first arrivals and neglects reflections which may have a significant impact on the illumination of the target.

In the computation of the wavenumber sampling, we rely on a single arrival high-frequency approximation. This means that multiple arrivals, such as reflections and multiples, are not considered in our quantification of the illumination at the target. The method would benefit from incorporating such effects, but extending the method to include multiple arrivals is far from being straightforward. This could be the scope of future work from this paper.

All things considered, we believe that our method fits monitoring applications, where a prior knowledge of the medium is assumed and the size of acquisition is usually small.

So far we have presented our method in light of the 2-D case, but there is no inherent difficulties that prevent its application to the 3-D case. The extension to 3-D and other developments will be the subject of a companion paper, and this represents a next step in the use of the method to realistic case scenarios, where 3-D acquisitions are today the standard in seismic design.

## DATA AVAILABILITY

No seismic data has been used, and the source code is available upon request to the authors.

## REFERENCES

* Amelin, N. & Yakovlev, A., 2024. Field trial results and potential application of autonomous marine drones, in *85th EAGE Annual Conference & Exhibition (Including the Workshop Programme)*, vol. 2024, pp. 1–5, European Association of Geoscientists & Engineers, doi: 10.3997/2214-4609.2024101460.

* Baraniuk, R.G. & Steeghs, P., 2017. Compressive sensing: a new approach to seismic data acquisition, *Leading Edge*, 36(8), 642–645.

* Barber, C.B., Dobkin, D.P. & Huhdanpaa, H., 1996. The quickhull algorithm for convex hulls, *ACM Trans. Math. Softw.*, 22(4), 469–483.

* Baudo, S., Calcagni, D., Esposito, A., Bettarelli, F., Ciccarone, T., Nini, C. & Hager, E., 2023. Acquiring hybrid High-resolution Streamer/OBN seismic for CCS in the complex environment of the Liverpool Bay, in

* 84th EAGE Annual Conference & Exhibition, vol. 2023, pp. 1–5, European Association of Geoscientists & Engineers. doi: 10.3997/2214-4609.202310902.

* Berger, J., Dutykh, D. & Mendes, N., 2017. On the optimal experiment design for heat and moisture parameter estimation, *Exp. Thermal Fluid Sci.*, 81, 109–122.

* Biondi, B.L., 2006. *3D Seismic Imaging*, Society of Exploration Geophysics.

* Blakely, R.J., 1995. *Potential Theory in Gravity and Magnetic Applications*, Cambridge Univ. Press,

* Chen, Y., Goetsch, P., Hoque, M.A., Lu, J. & Tarkoma, S., 2022. Dd-simplexed: adaptive delaunay triangulation for performance modeling and prediction on big data analytics, *IEEE Trans. Big Data*, 8(2), 458–469.

* Choksi, R. & Lu, X.Y., 2020. Bounds on the geometric complexity of optimal centroidal Voronoi tesselations in 3D, *Commun. Math. Phys.*, 377(3), 2429–2450.

* Curtis, A., 1999. Optimal experiment design: cross-borehole tomographic examples, *Geophys. J. Int.*, 136(3), 637–650.

* Curtis, A. & Snieder, R., 1997. Reconditioning inverse problems using the genetic algorithm and revised parametrization, *Geophysics*, 62(5), 1524–1532.

* Du, Q., Faber, V. & Gunzburger, M., 1999. Centroidal voronoi tessellations: applications and algorithms, *SIAM Rev.*, 41(4), 637–676.

* Fichtner, A. & Hofstede, C., 2022. A simple algorithm for optimal design in distributed fibre-optic sensing, *Geophys. J. Int.*, 233(1), 229–233.

* Gillies, S., van der Wel, C., Van den Bossche, J., Taves, M.W., Arnott, J. & Ward, B.C. & others, 2024. Shapely.

* Hardt, M. & Scherbaum, F., 1994. The design of optimum networks for aftershock recordings, *Geophys. J. Int.*, 117(3), 716–726.

* Harris, C.R. et al., 2020. Array programming with NumPy, *Nature*, 585(7825), 357–362.

* Karasözen, E. & Karasözen, B., 2020. Earthquake location methods, *GEM - Int. J. Geomath.*, 11(1), 13.

* Khinkis, L.A., Krzyzanski, W., Jusko, W.J. & Greco, W.R., 2009. D-optimal designs for parameter estimation for indirect pharmacodynamic response models, *J. Pharmacokinet. Pharmacodyn.*, 36(6), 523–539.

* Kijko, A., 1977. An algorithm for the optimum distribution of a regional seismic network-II. An analysis of the accuracy of location of local earthquakes depending on the number of seismic stations, *Pure appl. Geophys.*, 115(4), 1011–1021.

* Krampe, V., Edme, P. & Maurer, H., 2020. Optimized experimental design for seismic full waveform inversion: A computationally efficient method including a flexible implementation of acquisition costs, *Geophys. Prospect.*, 69(1), 152–166.

* Le Bouteiller, P., 2018. Eulerian Approach of Hamilton–Jacobi Equation with a Discontinuous Galerkin Method in Heterogeneous Anisotropic Medium: Application to Seismic Imaging, Ph.D. thesis, University Grenoble Alpes.

* Liner, C.L., Underwood, W.D. & Gobeli, R., 1999. 3-D seismic survey design as an optimization problem, *Leading Edge*, 18(9), 1054–1060.

* Liu, Y., Wang, W., Lévy, B., Sun, F., Yan, D.-M., Lu, L. & Yang, C., 2009. On centroidal Voronoi tessellation-energy smoothness and fast computation, *ACM Trans. Graph.*, 28(4), 101:1–101:17.

* Mancini, F. et al., 2023. A new generation of autonomous ocean-bottom nodes: their development and applications, in *84th EAGE Annual Conference & Exhibition*, vol. 2023, pp. 1–5, European Association of Geoscientists & Engineers. doi:10.3997/2214-4609.202310826.

* Martin, G.S., Wiley, R. & Marfurt, K.J., 2006. Marmousi2: an elastic upgrade for Marmousi, *Leading Edge*, 25(2), 156–166.

* Métivier, L. & Brossier, R., 2016. The seiscope optimization toolbox: a large-scale nonlinear optimization library based on reverse communication, *Geophysics*, 81(2), F11–F25.

* Mirebeau, J.-M. & Portegies, J., 2018. Hamiltonian fast marching: A numerical solver for anisotropic and non-holonomic eikonal PDEs, *ACM Trans. Graph*, 28, doi: 10.5201/ipol.2019.227.

* Mitchell, T.J., 2000. An algorithm for the construction of “D-Optimal” experimental designs, *Technometrics*, 42(1), 48–54.
Mora, P.R., 1989. Inversion = migration + tomography, *Geophysics*, **54**(12), 1575-1586.

Morrice, D.J., Kenyon, A.S. & Beckett, C.J., 2001. Optimizing operations in 3-D land seismic surveys, *Geophysics*, **66**(6), 1818-1826.

Nocedal, J. & Wright, S.J., 2006. *Numerical Optimization*, Springer, 2nd edn.

Ourabah, A., Popham, M., Einchcomb, C. & Goujon, N., 2020. A brief history of a disruptive nodal technology, in *1st Geoscience & Engineering in Energy Transition Conference*, vol. **2020**, pp. 1-5, European Association of Geoscientists & Engineers. doi: 10.3997/2214-4609.202021075.

Pagliccia, B., Dalton, K., Walker, C., Elder, K. & Jenneskens, R., 2018. METIS hits the ground in papua New Guinea, a field-proof innovative method to revolutionize onshore seismic acquisition, in *80th EAGE Conference and Exhibition 2018*, vol. **2018**, pp. 1-5, European Association of Geoscientists & Engineers. doi:

Pronzato, L., 2008. Optimal experimental design and some related control problems, *Automatica*, **44**(2), 303-325.

Rabinowitz, N. & Steinberg, D.M., 1990. Optimal configuration of a seismic network: a statistical approach, *Bull. seism. Soc. Am.*, **80**(1), 187-196.

Ridyard, D., Hite, D., Monk, D. & Mattox, B., 2023. A novel approach to OBN survey design for full waveform inversion, in *84th EAGE Annual Conference & Exhibition, Vol. 2023*, (1), pp. 1-5, European Association of Geoscientists and Engineers. doi:10.3997/2214-4609.202310190.

Sirgue, L. & Pratt, R.G., 2004. Efficient waveform inversion and imaging: a strategy for selecting temporal frequencies, *Geophysics*, **69**(1), 231-248.

Smith, K., 1918. On the standard deviations of adjusted and interpolated values of an observed polynomial function and its constants and the guidance they give towards a proper choice of the distribution of observations, *Biometrika*, **12**(1/2), 1-85.

Smucker, B., Krzywinski, M. & Altman, N., 2018. Optimal experimental design, *Nat. Methods*, **15**(8), 559-560.

Smythe, D.K., Smithson, S.B., Humphreys, C., Gillen, C., Kristoffersen, Y., Karaev, N.A., Garipov, V.Z. & Pavlenkova, N.I., 1994. Project images crust, collects seismic data in world's largest borehole, *EOS, Trans. Am. geophys. Un.*, **75**(41), 473-476.

Stork, C., 2019. *Global Land Seismic Acquisition Optimization by Accounting for Varying Noise, Obstacles, Non-Uniform Placement Costs, and Signal*, pp. 152-156. SEG Technical Program Expanded Abstracts 2019, Society of Exploration Geophysicists.

Stork, C., 2021. *Optimizing Seismic Land Acquisition and Processing for Surface Scattering Noise*, pp. 105-109. First International Meeting for Applied Geoscience \& Energy Expanded Abstracts, Society of Exploration Geophysicists.

Virieux, J. & Operto, S., 2009. An overview of full waveform inversion in exploration geophysics, *Geophysics*, **74**(6), WCC1-WCC26.

Virieux, J., Asnaashari, A., Brossier, R., Métivier, L., Ribodetti, A. & Zhou, W., 2017. An introduction to full waveform inversion, in *Encyclopedia of Exploration Geophysics*, pp. R1-1-R1-40, eds, Grechka, V. & Wapenaar, K., Society of Exploration Geophysics.

Virtanen, P. et al., 2020. SciPy 1.0: fundamental algorithms for scientific computing in python, *Nat. Methods*, **17**, 261-272.

Wang, S., Song, P., Tan, J., Xia, D., Du, G. & Wang, Q., 2024. Deblending of seismic data in the wavelet domain via a convolutional neural network based on data augmentation, *Geophys. Prospect.*, **72**(1), 213-228.

Wu, R.S. & Toksöz, M.N., 1987. Diffraction tomography and multisource holography applied to seismic imaging, *Geophysics*, **52**, 11-25.

Yang, P., Brossier, R., Métivier, L., Virieux, J. & Zhou, W., 2018. A time-domain preconditioned truncated Newton approach to multiparameter visco-acoustic full waveform inversion, *SIAM J. Sci. Comput.*, **40**(4), B1101-B1130.