System and method of simulating seismic wave propagation

Abstract

A method of simulating seismic wave propagation in media, the method implemented by a computer comprising a processor and a memory, the method comprising: modelling, using an earth model having a grid, a medium as a set of elastic tensors, wherein respective elastic tensors of the set thereof have tilted orientations relative to the grid; and simulating seismic wave propagation in the medium using a wave-equation operator, wherein the wave-equation operator is invariant to orientation changes of the respective elastic tensors of the set thereof.

Description

SYSTEM AND METHOD OF SIMULATING SEISMIC WAVE PROPAGATIONCROSS-REFERENCE TO RELATED APPLICATIONS

[0001] The present application claims priority to U. S. Patent Application 19 / 079,968 filed on March 14, 2025, which is incorporated by reference in its entirety herein.FIELD

[0002] The present disclosure relates to systems and methods of simulating seismic wave propagation, for example, for depth-imaging and inversion of earth model parameters using seismic data in the context of oil and gas (hydrocarbon) exploration.BACKGROUND

[0003] As hydrocarbon exploration and production move into more remote, deeper, and challenging geological settings, there is a growing demand for better and more reliable numerical algorithms to produce the necessary images and earth models.

[0004] The trend in both depth-imaging and inversion of seismic data over the past few decades has been the inclusion of more physical wave phenomena to match the growing complexity of the earth models. Among these is the increased use of anisotropy of the physical medium in which the seismic wave field propagates to improve the match between real input data and data simulations. Currently, the most cutting-edge numerical simulations of seismic wave propagation assume that the underlying earth models are composed of sedimentary rocks that can be described by elastic tensors with Tilted Transverse Isotropy (TTI) or Tilted Orthorhombic (TORT) symmetries, both often observed throughout the geologic sections of hydrocarbon-bearing basins and plays. Despite the many advances over the last 20 years in the accuracy of simulating seismic wave, the wave-equation operators used in these simulations still suffer from numerical instabilities. Addressing this issue is crucial for the success of seismic imaging and inversion in complex geological settings and thus to the success of hydrocarbon exploration and production.

[0005] Hence, there is a need to improve methods of simulating seismic wave propagation.SUMMARY

[0006] A first aspect provides a method of simulating seismic wave propagation in media, the method implemented by a computer comprising a processor and a memory, the method comprising:

[0007] modelling, using an earth model having a grid, a medium as a set of elastic tensors, wherein respective elastic tensors of the set thereof have tilted orientations relative to the grid; and

[0008] simulating seismic wave propagation in the medium using a wave-equation operator, wherein the wave-equation operator is invariant to orientation changes of the respective elastic tensors of the set thereof.

[0009] A second aspect provides a method of simulating wave propagation in media, the method implemented by a computer comprising a processor and a memory, the method comprising:

[0010] modelling, using a model having a grid, a medium as a set of elastic tensors, wherein respective elastic tensors of the set thereof have tilted orientations relative to the grid; and

[0011] simulating wave propagation in the medium using a wave-equation operator, wherein the wave-equation operator is invariant to orientation changes of the respective elastic tensors of the set thereof.

[0012] A third aspect provides a computer, comprising a processor and a memory, configured to perform the method according to any of the first aspect and / or the second aspect;

[0013] a computer program comprising instructions which, when executed by a computer, comprising a processor and a memory, cause the computer to perform the method according to any of the first aspect and / or the second aspect; and / or

[0014] a tangible (optionally, a non-transitory) computer-readable recording medium storing instructions which, when executed by a computer, comprising a processor and a memory, cause the computer to perform the method according to any of the first aspect and / or the second aspect.BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Fora better understanding of the disclosure, and to show how implementations of the same may be brought into effect, reference will be made, by way of example only, to the accompanying diagrammatic Figures, in which:

[0016] Figure 1 illustrates a cylinder and a Moebius strip.

[0017] Figure 2 illustrates a frame bundle concept forTTI and TORT models.

[0018] Figure 3 illustrates a method for simulating.

[0019] Figure 4 illustrates an example computing system that may implement various aspects of the disclosure.

[0020] Figure 5 illustrates a numerical experiment on a model with transverse isotropic symmetrical axis changing from vertical to horizontal.

[0021] Figure 5 to 12 illustrates a simulation example using a simple 2D velocity model, demonstrating that the orientation-invariant wave equations described herein can match or outperform traditional methods.DETAILED DESCRIPTIONSimulating seismic wave propagation

[0022] The first aspect provides a method of simulating seismic wave propagation in media, the method implemented by a computer comprising a processor and a memory, the method comprising:

[0023] modelling, using an earth model having a grid, a medium as a set of elastic tensors, wherein respective elastic tensors of the set thereof have tilted orientations relative to the grid; and

[0024] simulating seismic wave propagation in the medium using a wave-equation operator, wherein the wave-equation operator is invariant to orientation changes of the respective elastic tensors of the set thereof.

[0025] In this way, the simulation is improved compared with conventional simulations because the wave-equation operator is invariant to orientation changes of the respective elastic tensors of the set thereof, for example by accounting for and / or including respective orientations of the elastic tensors of the set thereof. In this way, the simulation is improved compared with conventional simulation, particularly where there are rapid and / or short-range spatial variations in tilt of the respective elastic tensors of the set thereof, for example in earth models having complex geological settings such as near reservoir bounding faults or salt bodies. In contrast, conventional simulations of such models typically suffer from numerical instabilities. In this way, imaging of earth models (e.g., Reverse Time Migration, RTM) and / or inference of earth model properties (e.g., Full Waveform Inversion, FWI, or in least-squares Reverse-time Migration, LSRTM) may be improved.Wave-equation operators

[0026] Conventional wave-equation operators designed for media having Tilted Transverse Isotropy (TTI) and Tilted Orthorhombic (TORT) symmetries have become standard in conventional seismic imaging for complex geological settings. These conventional wave-equation operators have been widely used in seismic exploration since the early 2000s, following the introduction of pseudoacoustic approximations. Extensive research and development, particularly in the 2010s, have made these conventional wave-equation operators integral to Reverse Time Migration (RTM) and Full-Waveform Inversion (FWI) algorithms in both academic and industrial applications.

[0027] A common feature of these conventional wave-equation operators is the use of directional derivatives to account for the spatial variation of TTI and TORT symmetry axes. Despite their widespread application and success, directional derivatives have limitations because they are not truly covariant differential operators in general. Consequently, numerical instabilities arise, particularly when dealing with pronounced spatial variations in the tilt of symmetry axes. Various ad-hoc solutions, such as smoothing the tilt fields before simulations, have been proposed, but these approaches are not always satisfactory, especially in complex geological settings where rapid spatial variations in tilt are common and need to be considered (e.g., near reservoir bounding faults or salt bodies). These solutions often fail to address the root of the problem, leading to persistent challenges in finding stable adjoint operators for TTI and TORT media.

[0028] The method according to the first aspect addresses and / or overcomes the limitations of the conventional use of directional derivatives by providing a new mathematical framework that inherently accounts for the symmetries of the elastic tensor, for example by adding tensor orientation as another set of coordinates to the problem of seismic wave propagation.

[0029] Particularly, the method according to the first aspect may: transform velocity models into curved manifolds; provide adjoint differential operators for propagating waves in such curved manifolds; and / or recognize that the variable tilt of symmetry elements of the elastic tensor has physical meaning and is a direct manifestation of the curvature of the modelling space itself. This new framework leads to stable wave-equation operators for TTI and TORT symmetries, for example.

[0030] The disclosure provides a novel representation of heterogeneous and anisotropic earth models using mathematical structures called frame bundles, for example. In this representation, local frames align with the symmetry elements (axes or planes) of the elastic tensor at each point in the earth model.

[0031] This approach allows for the incorporation and leverage of the inherent symmetries of the elastic tensors characterizing an earth model. Consequently, wave-equation operators can be written that are invariant to spatial changes in the orientation of symmetry elements of the elastic tensor, which may have Tilted Transverse Isotropy (TTI), Tilted Orthorhombic (TORT), or Tilted Monoclinic symmetries, for example.

[0032] The disclosure provides simplified wave equations but endows their differential operators with a key attribute: it makes them compatible with the metric of the resulting frame-bundle manifold representing earth models. This compatibility ensures true adjoint operators even in the presence of rapid spatial variation of symmetry axes. This is a significant improvement over current state-of-the-art wave-equation operators, leading to more accurate and reliable seismic imaging and inversion results, as current operators for TTI and TORT media lack this property.

[0033] The method according to the first aspect offers a more intuitive and physically meaningful way to parameterize anisotropic earth models compared with conventional methods. The method according to the first aspect provides new physical insights into seismic wave phenomena in anisotropic and heterogeneous media because it connects topology to elastic tensor symmetry orientation. The earth models produced by the method according to the first aspect are mathematically coherent manifolds which are in general curved, with curvature (a topological property) proportional to the spatial variation of symmetry orientations. This topological information should be valuable for designing better inversion algorithms, because it provides unprecedented constraints on the estimations of symmetry orientations, which are typically assumed to be orthogonal to reflectors on depth-migrated seismic images, for example.

[0034] The method according to the first aspect is applicable to seismic imaging and inversion of seismic data in various settings, from exploration and production of hydrocarbons to monitoring of subsurface processes like CO2 sequestration or hydraulic fracturing. The method according to the first aspect is general enough to be applied in other fields of science and engineering where wave phenomena in anisotropic elastic media are of interest.

[0035] The method according to the first aspect is also flexible to be used with various numerical methods, such as finite-difference, finite-element, and spectral-element methods, for solving waveequation partial differential equations in pseudo-acoustic or elastic media. The method according to the first aspect may be employed for imaging methods (e.g., Reverse Time Migration, RTM) and / or inference of earth model properties (e.g., Full Waveform Inversion, FWI, or in least-squares Reversetime Migration, LSRTM).Method

[0036] The method includes simulating seismic wave propagation in media. Generally, a solid body may be deformed by the application of an external force. If the solid body is perfectly elastic, the solid body will return to its original shape once that external force is removed. In the context of exploration seismology, the earth can generally be considered as perfectly elastic because the stresses generated by seismic exploration activities are too small to permanently deform subsurface rocks. When an impulsive or transitory stress is applied to a finite area on the surface of an elastic solid body, a strain is generated in the immediately adjacent sub-volume. The strained sub-volume then transfers stress in turn to mutually adjacent sub-volumes within the solid body, which generates strains in these mutually adjacent sub-volumes. In this way, an impulsive stress propagates through a solid body as an elastic wave. Elastic waves that propagate in the earth are known as seismic waves. Generally, a seismic wave is a mechanical wave of elastic energy that travels through the earth or another planetary body. The propagation velocity of a seismic wave depends on the density and the elasticity of the solid body as well as the type of wave. In one example, the media comprise and / or are solid bodies, for example planets, such as the earth, or parts, such as volumes, thereof. In one example, the media comprise and / or are elastic media, for example perfectly elastic media and / or considered to be perfectly elastic media. In one example, the media comprise and / or are elastic solid bodies, for example perfectly elastic solid bodies and / or considered to be perfectly elastic solid bodies, for example planets, such as the earth, or parts, such as volumes, thereof. In one example, the media comprise and / or are heterogeneous and / or anisotropic media, for example as described with respect to the medium mutatis mutandis.Computer

[0037] The method is implemented by the computer comprising the processor and the memory. In other words, the method comprises and / or is a computer-implemented method.Modelling

[0038] The method comprises modelling, using the earth model having the grid, the medium as the set of elastic tensors (also known as elasticity tensors, elastic modulus tensors and stiffness tensors). Modelling, using an earth model having a grid, the medium as a set of elastic tensors is known, for example using a 2D or a 3D grid having points or nodes. Generally, an elasticity tensor is a fourthrank tensor describing the stress-strain relation in a linear elastic material. A grid is also known as acomputational grid or a coordinate system. The grid is used to represent, at least in part, the earth model and / or a part thereof.

[0039] The respective elastic tensors of the set thereof have tilted orientations relative to the grid. In this way, symmetry elements of different tensors at each point or node of the grid may be oriented differently, for example mutually differently. In one example, two or more mutually adjacent elastic tensors of the set thereof, for example at mutually adjacent points or nodes of the grid, have mutually different tilt orientations. In one example, for TTI symmetry, a symmetry element comprises an axis, for example a single axis, of rotation invariance by any angle (or equivalently, a plane of isotropy). In one example, for TORT media, symmetry elements comprise three mutually orthogonal axes of rotation invariance to angles equal to 180 degrees (or equivalently to three mutually perpendicular mirror symmetry planes).

[0040] In one implementation, the symmetry elements are spatially variable, for example varying from pointto point or node to node on a grid such as a spatial coordinate grid, such as having different directions or orientations (i.e. tilts) and / or different magnitudes from point to point or node to node on a grid such as a spatial coordinate grid.Earth model

[0041] Generally, an earth model is a representation of the Earth's internal structure, typically including parameters such as density, seismic P-waves velocity (a), and S-waves velocities (P).

[0042] In one implementation, the earth model comprises and / or is a heterogeneous earth model, for example comprising one or more heterogeneities having mutually different elasticities and / or densities, whereby respective propagation velocities of a seismic wave therein and / or therethrough are mutually different, such as one or more heterogeneous layers, one or more reservoirs such as one or more hydrocarbon reservoirs, one or more faults such as reservoir bounding faults and / or one or more included bodies such as salt bodies.

[0043] In one implementation, the earth model comprises and / or is an anisotropic earth model, for example comprising directionally dependent physical properties such as elasticity, whereby respective propagation velocities of a seismic wave are directionally dependent. Generally, anisotropy stems from the effects of inherent properties of rocks (e.g., like the alignment of phyllosilicate minerals) or from the combined effect of various thin layers of isotropic rocks in a sedimentary formation on the propagation of seismic wavefields. For example TTI symmetries areusually associated with thick sequences of shales in sedimentary basis, which show tilt in relation to the earth’s surface because of various episodes of faulting or folding throughout the geologic time.

[0044] In one implementation, the earth model comprises and / or is a heterogeneous, anisotropic earth model.Medium

[0045] In one implementation, the medium comprises and / or is a solid body, for example a planet, such as the earth, or a part, such as a volume, thereof.

[0046] In one implementation, the medium comprises and / or is an elastic medium, for example a perfectly elastic medium and / or considered to be a perfectly elastic medium.

[0047] In one implementation, the medium comprises and / or is an elastic solid body, for example a perfectly elastic solid body and / or considered to be a perfectly elastic solid body, for example a planet, such as the earth, or a part, such as a volume, thereof.

[0048] In one implementation, the medium comprises and / or is a heterogeneous medium, for example comprising one or more heterogeneities having mutually different elasticities and / or densities, such as and / or due to mutually different textures, mineral composition, permeabilities and / or porosities, whereby respective propagation velocities of a seismic wave therein and / or therethrough are mutually different, such as one or more heterogeneous layers, one or more reservoirs such as one or more hydrocarbon reservoirs, one or more faults such as reservoir bounding faults and / or one or more included bodies such as salt bodies.

[0049] In one implementation, the medium comprises and / or is an anisotropic medium, for example comprising directionally dependent physical properties such as elasticity, such as whereby respective propagation velocities of a seismic wave are directionally dependent.

[0050] In one implementation, the medium comprises and / or is a heterogeneous, anisotropic medium.

[0051] In one implementation, the medium has Tilted Transverse Isotropy, TTI, symmetry, Tilted Orthorhombic, TORT, symmetry and / or Tilted Monoclinic symmetry, preferably TTI and / or TORT. Other symmetries are known.Simulating

[0052] In an implementation, the method comprises simulating seismic wave propagation in the medium using the wave-equation operator, wherein the wave-equation operator is invariant to orientation changes of the respective elastic tensors of the set thereof.

[0053] In this way, the simulation is improved compared with conventional simulations because the wave-equation operator is invariant to orientation changes of the respective elastic tensors of the set thereof, for example by accounting for and / or including respective orientations of the elastic tensors of the set thereof. In this way, the simulation is improved compared with conventional simulation, particularly where there are rapid and / or short-range spatial variations in tilt of the symmetry elements, for example in earth models having complex geological settings such as near reservoir bounding faults or salt bodies. In contrast, conventional simulations of such models typically suffer from numerical instabilities.

[0054] In one implementation, simulating seismic wave propagation in the medium using the waveequation operator comprises aligning respective local frames of the wave-equation operator with the respective symmetry elements of the respective elastic tensors of the set thereof, for example at each point or node on a grid such as a spatial coordinate grid. In this way, the wave-equation operator is invariant to orientation changes of the respective elastic tensors of the set thereof.Fiber bundles

[0055] In one implementation, aligning respective local frames of the wave-equation operator with the respective symmetry elements comprises providing the respective local frames by a fiber bundle. The fiber bundles provide a framework for anisotropic velocity models, for example, by generalizing Cartesian products, allowing the creation of complex manifolds, such as Riemannian manifolds describing elastic tensors by three spatial coordinates and three local orientation coordinates, by combining simpler ones, as described herein.

[0056] As detailed below, the fiber bundle approach splits the tangent space of the total space into two complementary pieces: one piece along the fibers and another piece perpendicular (more generally: complementary) to the fibers. Consequently, the connection operator (also known as the covariant derivative) acting on a vector or tensor field T defined over the total space manifold and taken along a vector X is given by the summation of two operators: DXT = VxT + E u(X)T, where Vrepresents the part of the connection accounting for changes in position in the base space, and w for changes in the orientation of the elastic tensor (i.e., in the local frame orientations). The connection 1-form, M, acts over each index of T, thus the summation symbol in the definition of DXT. It should be understood that (differential) 1-forms are dual to vectors, i.e., on manifolds, 1-forms represent linear functions acting on vectors to produce scalars. For simplicity and to highlight the contribution of spatial variation of orientations of anisotropic axes into wave propagation, the base space is taken to represent a Euclidean space, such that Vx= dx, i.e., changes along base space directions are given by directional derivatives. Technical details aside, the operator M in the frame bundle has a very specific form: it is a skew-symmetric 3 x 3 matrix, whose elements are given by < > = where h is the rotation / attitude matrix determining the orientation of the elastic tensors at each point in the anisotropic velocity model. Thus, M is a model parameter, independent of propagation, which can be precomputed or obtained on the fly, once the attitude matrix is set, regardless of its parameterization using Euler angles or quaternions.

[0057] Given these preliminary theoretical considerations, any wave equation operator for media may be recast using the proposed frame-bundle approach. This may involve replacing any partial or directional derivatives with the bundle connection D and ensuring the resulting equation is frame independent. For example, starting from the wave equation (1):d j = A'‘ dtdrakr

[0058] (1)

[0059] may be recast as equation (2):d j = AlfjD Drffkr= A^D(tr(Da))fc;= A (tr^^ )))^ I0060] (2)

[0061] where A is the density-normalized elastic tensor and repetition of an index in one upstairs and one downstairs position indicates a summation over it (i.e., using Einstein’s summation convention). After the second equality sign in equation 2, we switched to a more abstract notation to show that the action of the two connection operators D results in the gradient of the divergence of the stress tensor a. Because D is compatible with the metric, it commutes with the trace operator and so one reaches the expression with the term D(Do-). This is not a tensorial quantity, thus not frame independent. It is, nevertheless, part of the Hessian of the stress tensor, £>y<7 = Dx(DYa) -DDxYa, which is tensorial thanks to the extra term DDXYU (with X and Y being vectors in the manifold). Because equation (1), the starting point, is set up in a Euclidean manifold (where the second term of the Hessian is identically zero) the translation from directional derivatives to covariant ones requires that dldsakrbe understood as the full Hessian. Thus, equation (2) may be written as equation (3):= A^tr(D2a)]lcl= -A^alcl

[0062] (3)

[0063] where A<rklrepresents the Laplacian of the stress tensor, which, however, is not simply the summation of second partial derivatives. In curved spaces, like the frame bundle, second covariant derivatives do not commute, but their difference defines a crucial quantity that measures the curvature of the total space itself, the curvature tensor R X, Y), i.e. equation (4):£)X, Y(T — =DXDY< T — DYDxa — D[XY](T = R(X, Y)a

[0064] (4)

[0065] with [X, y] = dxY - dYX also a vector. Equation (4) underscores that the Laplacian in equation (3) contains terms related to the curvature tensor, since explicitly from equation (5):A<r = — tr [DXY< T] = — tr[Z)yX<r + R (X, y)cr]

[0066] (5)

[0067] This is a physically insightful result: this equation shows that space curvature leads to extra accelerations in wave-equations because of elastic tensor misalignments. Here, curvature is created by misalignment of elastic tensors (in general relativity, by the presence of mass). Also, the wave equation describes how stress-waves propagate (the acceleration term) as a function of spatial changes in stress distribution (the Laplacian term). And the equation shows how the misalignment of elastic tensors creates extra terms in these spatial changes in stress, captured by the Laplacian. In more detail, firstly, it clarifies why the directional derivative approach runs into problems, because the resulting Laplacian operators would then not have terms with the connection 1-form < D (found in D) and the curvature R. Secondly, it links differences in orientation of symmetry axes in models, for example, to extra accelerations in the wave equation, coming from the action of the curvature (tensor) on the stress perturbations describing the waves. The following expression for the Laplacian of second-order tensors in Riemannian manifolds, is given by equation (6):(Acr)fci= + Richer / + Ric^ - 2Rijkla^,

[0068] (6)

[0069] where Ric is the Ricci tensor, obtained from contraction of two indices of the curvature tensor R. The Lichnerowicz Laplacian in equation (6) is self-adjoint and it neatly separates the curvature effects when computing Laplacians. The term DlD.a!, the rough Laplacian, is what one would get if the manifold were flat. It comprises the familiar d’ terms plus terms with the connection 1-form co. Note that the Laplacian returns a tensor of the same dimension and having the same symmetriesof the input tensor. In the frame bundle, this means that A<r can also be represented by a 3 x 3 symmetric matrix. Both the curvature and the Ricci tensors are evaluated to where the wavefield is, making finite-difference stencils redundant and thus having negligible impact on numerical computations during propagation. The curvature tensor R can be calculated from the metric tensor of the bundle but is more efficiently computed from the connection 1-form. Indeed, the covariant derivative of > results in another skew-symmetric 3 x 3 matrix, given by equation (7):n(X, T) = D X. Y) = dxM(Y) - dYa)(X) + to(X) A(o(7)

[0070] (7)

[0071] A standard result of differential geometry then connects the matrix (1 to the curvature tensor R, because each element of this matrix Q is a function of the components of the curvature tensor R in a local coordinate frame xl, given by equation (8):

[0073] Both equations (7) and (8) include the wedge product A, the product operation for differential forms. The wedge product A is the tensor product followed by an antisymmetrization operation, which consists of summation over all components of the tensor product, but applying alternating signs based on the permutation of the indices in each term. The antisymmetrization ensures that the resulting differential / c-forms are always alternating functions of k input vectors, thus encoding in them oriented lines, planes, volumes, and their higher-dimensional counterparts, enabling and streamlining calculus (and thereby wave propagation) in general manifolds. The wedge product A allows for the combination of lower-dimensional forms into higher-dimensional ones. For example, wedging two 1 -forms produces a 2-form, an alternating bilinear function that takes two vectors and outputs the oriented area spanned by them in the plane of the 2-form, with the sign of the result changing whenever one swaps the order of the input vectors. Thus, in equation (8), the elements of the matrix n are summations of 2-forms whose coefficients are the components of the curvature tensor R and the unit areas dxkAdx1are defined by wedging the 1 -forms dxkand dxlrelated to a set of coordinates xl(i = 1,used to locally parametrize the frame bundle.

[0074] With all the elements above, one is now able to write covariant wave equations for media, as above, going from equation (1) to equation (3), for example. The key point is that defining a frame bundle parameterized by the local orientation of elastic stiffness tensors provides wave equations that are invariant to changes in the orientation of the symmetry axes of these same elastic tensors.

[0075] In one implementation, providing the respective local frames by the fiber bundle comprises attaching a copy of an orientation grid (for example, of the respective symmetry elements) to respective points (or nodes) in a spatial coordinate grid of the earth model. In this way, elastic tensor orientation may be decoupled from the spatial geometry.

[0076] In one implementation, providing the respective local frames by the fiber bundle comprises representing space of all orientations of respective local frames (for example, describing the local coordinate axes) using the respective fibers. In this way, an abstraction is provided to enable deployment of the fiber bundles to understand elastic wave propagation in the model. In this way, a connection operator may be defined, as described herein, that is compatible with a metric of a 6D manifold, for example. This has two significant consequences: firstly, it ensures that one can obtain true adjoint operators for wave propagation in the frame bundle manifold; and secondly, it provides differential operators that automatically account for spatial changes in orientation of symmetry axes, for example TTI or TORT symmetry axes, when taking (covariant) derivatives, thereby simplifying equations and numerical implementations, as there is no longer need to do and track rotations of elastic tensors or of directional derivatives, for example.

[0077] In one implementation, providing the respective local frames by the fiber bundle comprises representing the respective local frames axially and complementary to the respective fibers. In this way, the fiber bundle splits the tangent space of the total space into two complementary pieces: one piece along the fibers and another piece perpendicular (more generally: complementary) to the fibers. Consequently, the connection operator (also known as the covariant derivative) acting on a vector or tensor field defined over the total space manifold and taken along a vector is given by the summation of two operators, as described herein.

[0078] In one implementation, providing the respective local frames by the fiber bundle comprises providing the respective local frames by the fiber bundle having sixdimensions, for example generally curved manifolds; wherein a base space provides coordinates for location in 3D space; and wherein respective fibers represent orientation manifolds, for example parameterized by coordinates and describing respective local orientations of the respective symmetry elements. In this way, by defining a frame bundle parameterized by the local orientation of elastic stiffness tensors, wave equations are provided that are invariant to changes in the orientation of the symmetry axes of these same elastic tensors.

[0079] Outputting, imaging, inference & adjusting

[0080] In one implementation, the method comprises outputting a result of simulating seismic wave propagation in the medium. In this way, the result of the simulating seismic wave propagation in the medium may be compared, for example, with a result of measuring seismic wave propagation in a physical medium (i.e. the medium that has been modelled, for example represented with a subsurface model).

[0081] In one implementation, the method comprises comparing the result of simulating seismic wave propagation in the medium with a result of measuring seismic wave propagation in a physical medium. In one implementation, the method comprises determining one or more differences based on a result of the comparing. In one implementation, the method comprises updating the model based on the one or more differences. In one implementation, the method comprises iteratively updating the model, for example by repeating the outputting, the comparing, the determining and / or the updating. In this way, the model may be adjusted.

[0082] Provided is a method of seismic imaging of an earth model, comprising the method according to the first aspect. In this way, the seismic imaging may be compared, for example, with seismic imaging of a physical medium (i.e. the medium that has been modelled, for example represented with a subsurface model).

[0083] Provided is a method of seismic inference of an earth model, comprising the method according to the first aspect. In this way, properties of a physical medium (i.e. the medium that has been modelled, for example represented with a subsurface model) may be inferred from the seismic inference.

[0084] Provided is a method of adjusting an earth model, comprising the method according to the first aspect.

[0085] Provided is a method of inversion of seismic data, comprising the method according to the first aspect.

[0086] Provided is a method where, in one implementation, the method provides a method of and / or in one example the method comprises seismic wave propagation modelling in heterogeneous and anisotropic elastic media with TTI or TORT symmetry, using fiber bundle theory to derive waveequations that are invariant to orientation changes of symmetry elements (axes or planes) throughout an earth model.

[0087] Provided is a method where, in one implementation, the method provides a method of and / or in one example the method comprises defining heterogeneous anisotropic velocity models as frame bundles of six dimensions that are generally curved manifolds, in which the base space provides the coordinates for location in 3D space, while the fibers represent orientation manifolds (e.g., the group of rotations in 3D), parameterized by another three coordinates and describing the local orientation of symmetry elements of the elastic stiffness tensor.

[0088] Provided is a method where, in one implementation, the method provides a method of and / or in one implementation the method comprises defining a connection operatorthat is compatible with the metric of the resulting frame bundle manifold, ensuring by design true adjoint differential operators in wave propagation and accounting for spatial changes in the orientation of TTI or TORT symmetry axes.

[0089] Provided is a method where, in one implementation, the method provides a method of and / or in one implementation the method comprises simplifying wave equations and their numerical implementations by using connection operators that automatically account for spatial changes in elastic tensor orientations, eliminating the need to track relative rotations / orientations of these tensors or of differential operators.

[0090] Provided is a method where, in one implementation, the method provides a method of and / or in one implementation the method comprises improving stability and reducing numerical artifacts in wavefield simulations (in comparison to state-of-the art methods using directional derivatives) by using the afore-mentioned connection operators and corresponding 2-form connection curvature to accurately account for spatial variation of orientation of elastic tensors in a velocity model.

[0091] Provided is a method where, in one implementation, the method provides a method of and / or in one implementation the method comprises lending physical meaning to orientation changes of symmetry elements as manifold curvature that creates new tidal forces in wave equations, providing new insights into the behavior of wavefields propagating in heterogeneous and anisotropic elastic media.

[0092] Provided is a method where, in one implementation the method provides a method of and / or in one implementation the method comprises introducing new model parameters constraints in the forms of a curvature tensor (or equivalently the 2-form connection curvature) and its various contractions (Ricci tensor and scalar curvature) for inversion problems, providing new ways to estimate symmetry orientations from seismic data.

[0093] Provided is a method where, in one implementation the method provides a method of and / or in one implementation the method comprises enhancing reflectivity and velocity model estimation in inverse problems such as LSRTM or FWI by using the frame bundle approach to obtain superior fit between observed and modelled data and to separate the effects of curvature stemming from underlying space geometry from those caused by misalignments of elastic tensor orientations.

[0094] Provided is a method where, in one implementation, the method provides a method of and / or in one implementation the method comprises developing more accurate and efficient solvers that can adaptively account for frame rotations, leveraging the frame bundle approach to improve wave propagation modelling in anisotropic media.

[0095] Provided is a method where, in one implementation, the method provides a method of and / or in one implementation the method comprises extending wave propagation modelling to tilted Monoclinic symmetry cases by using elastic tensors with monoclinic symmetry.

[0096] Provided is a method where, in one implementation, the method provides a method of and / or in one implementation the method comprises encouraging further investigation into orientation parameterization, such as using Euler angles or quaternions, to improve the accuracy and efficiency of wave propagation modelling in anisotropic media.Simulating wave propagation

[0097] The second aspect provides a method of simulating wave propagation in media, the method implemented by a computer comprising a processor and a memory, the method comprising:

[0098] modelling, using a model having a grid, a medium as a set of elastic tensors, wherein respective elastic tensors of the set thereof have tilted orientations relative to the grid; and

[0099] simulating wave propagation in the medium using a wave-equation operator, wherein the wave-equation operator is invariant to orientation changes of the respective elastic tensors of the set thereof.

[0100] The method according to the second aspect may be as described and / or include any step as described with respect to the first aspect, mutatis mutandis.

[0101] In one implementation, the media comprise and / or are solid bodies, for example biological bodies such as anatomical bodies or structural bodies such as manufactured components orassemblies, or parts, such as volumes, thereof. In one implementation, the media comprise and / or are elastic media, for example perfectly elastic media and / or considered to be perfectly elastic media. In one implementation, the media comprise and / or are elastic solid bodies, for example perfectly elastic solid bodies and / or considered to be perfectly elastic solid bodies, for example biological bodies such as anatomical bodies or structural bodies such as manufactured components or assemblies, or parts, such as volumes, thereof. In one example, the media comprise and / or are heterogeneous and / or anisotropic media, for example as described with respect to the medium mutatis mutandis.

[0102] Provided is a method of imaging of a model, comprising the method according to the second aspect.

[0103] Provided is a method of inference of a model, comprising the method according to the second aspect.

[0104] Provided is a method of adjusting a model, comprising the method according to the second aspect.

[0105] Provided is a method of inversion of data, comprising the method according to the second aspect.Computer, computer program and / or CRM

[0106] The third aspect provides a computer, comprising a processor and a memory, configured to perform the method according to any of the first aspect and / or the second aspect;

[0107] a computer program comprising instructions which, when executed by a computer, comprising a processor and a memory, cause the computer to perform the method according to any of the first aspect and / or the second aspect; and / or

[0108] a tangible (optionally, a non-transitory) computer-readable recording medium storing instructions which, when executed by a computer, comprising a processor and a memory, cause the computer to perform the method according to any of the first aspect and / or the second aspect.

[0109] The disclosure applies concepts of differential geometry to the problem of seismic wave propagation in heterogeneous anisotropic elastic media. The anisotropic velocity models are conceptualized as Riemannian manifolds described by six coordinates: three spatial coordinates and another three for the local orientation of the elastic tensor. This parameterization is locally trivial butglobally complex, akin to how the Earth’s surface appears locally flat (i.e., Euclidean) but is globally curved. Thus, anisotropic velocity models are not six-dimensional Euclidean spaces and instead require a more sophisticated framework. Fiber bundles provide this framework by generalizing Cartesian products, allowing the creation of complex manifolds by combining simpler ones. In this case, this combination (known as a fiber bundle) results from attaching a copy of the orientation grid (a manifold on its own) to every point in the spatial coordinate grid (the base space, also known as the base manifold). These fibers can "twist" as they attach to the base space, creating intricate global structures, which are stitched together smoothly by transition functions. The fiber bundle is characterized by the combined total space, the base space, and a projection map that matches points from the total space to the base space. Classic illustrations of fiber bundles are a cylinder and a Moebius strip. Both the cylinder and the Moebius strip may be represented by the combination of a circle as a base space with a real line representing the fibers attached to each point of that circle. Locally, both the cylinder C and the Moebius strip M look like Euclidean spaces, but in the Moebius strip, the orientation of the fibers changes along or around the circle, resulting in a total space that is not a cylinder (see Figure 1). Referring to Figure 1, the cylinder and the Moebius strip are two fiber bundles created by attaching copies of the real line (the fibers) to a circle (the base space). Locally they can be represented by the same flat space, but their global topologies are widely different because how the fibers twist in the Moebius strip. Finally, connections extend the notion of directional derivatives to fiber bundles. In curved spaces, connections allow the parallel transport of vectors and tensors along curves in the base space, enabling their comparison and differentiation across different fibers.

[0110] In the description of TTI or TORT velocity models (more generally: velocity models) herein, the so-called frame bundle is used, where the fiber represents the space of all orientations of the frame describing the local coordinate axes. The total space model is characterized by a local frame at each point in the base space, taken to be the same frame in which elastic tensors are given by their smallest number of parameters (for example, five in the case of TTI and nine for TORT; see Figure 2). Referring to Figure 2, the 2-D surface, including a concavity V and a convexity X, represents the computational grid G, which determines the spatial discretization of the model using general curvilinear grids. The mutually perpendicular axes A1, A2, A3 show local frame orientations that are attached to the base space but align with the elastic tensor symmetry elements at each grid position, not to the geometry of the underlying base space. This way, elastic tensor orientation is decoupled from the underlying spatial geometry. This abstraction enables deployment of the powerful tools of fiber bundles to understand elastic wave propagation in anisotropic models in a novel way. At once, it allows for a connection operator to be defined that is compatible with the metric of the 6D manifold. This has two significant consequences: firstly, it ensures that one can obtain true adjointoperators for wave propagation in the frame bundle manifold; and secondly, it provides differential operators that automatically account for spatial changes in orientation of symmetry axes, for example TTI or TORT symmetry axes, when taking (covariant) derivatives, thereby simplifying equations and numerical implementations, as there is no longer need to do and track rotations of elastic tensors or of directional derivatives, for example.

[0111] The fiber bundle approach splits the tangent space ofthe total space into two complementary pieces: one piece along the fibers and another piece perpendicular (more generally: complementary) to the fibers. Consequently, the connection operator (also known as the covariant derivative) acting on a vector or tensor field T defined over the total space manifold and taken along a vector X is given by the summation of two operators: DXT = VXT + X wOOT, where V represents the part of the connection accounting for changes in position in the base space, and w for changes in the orientation ofthe elastic tensor (i.e., in the local frame orientations). The connection 1-form, acts over each index of T, thus the summation symbol in the definition of DXT. It should be understood that (differential) 1 -forms are dual to vectors, i.e., on manifolds, 1 -forms represent linear functions acting on vectors to produce scalars. For simplicity and to highlight the contribution of spatial variation of orientations of anisotropic axes into wave propagation, the base space is taken to represent a Euclidean space, such that Vx= dx, i.e., changes along base space directions are given by directional derivatives. Technical details aside, the operator a) in the frame bundle has a very specific form: it is a skew-symmetric 3 x 3 matrix, whose elements are given by = h~1dh, where h is the rotation / attitude matrix determining the orientation of the elastic tensors at each point in the anisotropic velocity model. Thus, M is a model parameter, independent of propagation, which can be precomputed or obtained on the fly, once the attitude matrix is set, regardless of its parameterization using Euler angles or quaternions.

[0112] Given the above preliminary theoretical considerations, any wave equation operator for media with TTI or TORT symmetries, for example, may be recast using the proposed frame-bundle approach. This may involve replacing any partial ordirectional derivatives with the bundle connection D and ensuring the resulting equation is frame independent. For example, starting from the wave equation written in terms of stress tensor a components and considering constant density, equation (1):dt^j = drffkr

[0113] (7) may be recast as equation (2):dt ^ij = A^DiDrakr= A^D(tr(Da))fc;= A (tr (p^Da)^

[0114] (8)

[0115] where A is the density-normalized elastic tensor and repetition of an index in one upstairs and one downstairs position indicates a summation over it (i.e., using Einstein’s summation convention). After the second equality sign in equation 2, we switched to a more abstract notation to show that the action of the two connection operators D results in the gradient of the divergence of the stress tensor a. Because D is compatible with the metric, it commutes with the trace operator and so one reaches the expression with the term D(Dcr). This is not a tensorial quantity, thus not frame independent. It is, nevertheless, part of the Hessian of the stress tensor, Dxo = Dx(pYo) -DDXYO, which is tensorial thanks to the extra term DDXY< T (with X and Y being vectors in the manifold). Because equation (1), the starting point, is set up in a Euclidean manifold (where the second term of the Hessian is identically zero) the translation from directional derivatives to covariant ones requires that dldsakrbe understood as the full Hessian. Thus, equation (2) may be written as equation (3):d^ij = A^[tr^D2a)]kl= -A^Aakl

[0116] (9)

[0117] where Aaklrepresents the Laplacian of the stress tensor, which, however, is not simply the summation of second partial derivatives. In curved spaces, like the frame bundle, second covariant derivatives do not commute, but their difference defines a crucial quantity that measures the curvature of the total space itself, the curvature tensor R (%, / ), i.e. equation (4):D^ycr— )r,x°'=DXDYU — DYDxa — D[X y]<7 = R(X, Y)ff

[0118] (10)

[0119] with [X, Y] - dxY - dYX also a vector. Equation (4) underscores that the Laplacian in equation (3) contains terms related to the curvature tensor, since explicitly from equation (5):ACT = — tr[D(y<r] = — tr[Zy xa + R(X, K)<r]

[0120] (11)

[0121] This is a physically insightful result: this equation shows that space curvature leads to extra accelerations in wave-equations because of elastic tensor misalignments. Like in General relativity, where Gravity is a force stemming from the curvature of space-time. Here, curvature is created by misalignment of elastic tensors (in general relativity, by the presence of mass). Also, the wave equation describes how stress-waves propagate (the acceleration term) as a function of spatial changes in stress distribution (the Laplacian term). And the equation shows how the misalignment of elastic tensors creates extra terms in these spatial changes in stress, captured by the Laplacian. Inmore detail, firstly, it clarifies why the directional derivative approach runs into problems, because the resulting Laplacian operators would then not have terms with the connection 1-form co (found in D) and the curvature 7?. Secondly, it links differences in orientation of symmetry axes in TTI and TORT models, for example, to extra accelerations (akin to tidal forces) in the wave equation, coming from the action of the curvature (tensor) on the stress perturbations describing the waves. The following expression for the Laplacian of second-order tensors in Riemannian manifolds, as given by equation (6):(A<r)fc / = -D'DiOfci + Richer + RiciZo - 2Rijkicjli,

[0122] (12)

[0123] where Ric is the Ricci tensor, obtained from contraction of two indices of the curvature tensor R. The Lichnerowicz Laplacian in equation (6) is self-adjoint and it neatly separates the curvature effects when computing Laplacians. The term DlD.a!, the rough Laplacian, is what one would get if the manifold were flat. It comprises the familiar d’terms plus terms with the connection 1-form co. Note that the Laplacian returns a tensor of the same dimension and having the same symmetries of the input tensor. In the frame bundle, this means that Ao- can also be represented by a 3 x 3 symmetric matrix. Both the curvature and the Ricci tensors are evaluated to where the wavefield is, making finite-difference stencils redundant and thus having negligible impact on numerical computations during propagation. The curvature tensor 7? can be calculated from the metric tensor of the bundle but is more efficiently computed from the connection 1-form. Indeed, the covariant derivative of co results in another skew-symmetric 3 x 3 matrix, given by equation (7):n X, y) = (Dw)(X, r) = dxcj(Y) - dY(o(X) + co(X) A co(Y)

[0124] (7)

[0125] A standard result of differential geometry then connects the matrix n to the curvature tensor R, because each element of this matrix Q is a function of the components of the curvature tensor R in a local coordinate frame xl, given by equation (8):

[0126] (8)nj = 2 ^jki dx>cA

[0127] Both equations (7) and (8) include the wedge product A, the product operation for differential forms. The wedge product is the tensor product followed by an antisymmetrization operation, which consists of summation over all components of the tensor product, but applying alternating signs based on the permutation of the indices in each term. The antisymmetrization ensures that the resulting differential / c-forms are always alternating functions of k input vectors, thus encoding inthem oriented lines, planes, volumes, and their higher-dimensional counterparts, enabling and streamlining calculus (and thereby wave propagation) in general manifolds. The wedge product allows for the combination of lower-dimensional forms into higher-dimensional ones. For example, wedging two 1 -forms produces a 2-form, an alternating bilinear function that takes two vectors and outputs the oriented area spanned by them in the plane of the 2-form, with the sign of the result changing whenever one swaps the order of the input vectors. Thus, in equation (8), the elements of the matrix (1 are summations of 2-forms whose coefficients are the components of the curvature tensor R and the unit areas dxk / \dxlare defined by wedging the 1 -forms dxkand dxlrelated to a set of coordinates (i = 1,used to locally parametrize the frame bundle.

[0128] With all the elements above, one is now able to write covariant wave equations for TTI and TORT media, for example, as above, going from equation (1) to equation (3), for example. The key point is that defining a frame bundle parameterized by the local orientation of elastic stiffness tensors provides wave equations that are invariant to changes in the orientation of the symmetry axes of these same elastic tensors.

[0129] The disclosure herein provides a new, mathematically rigorous method to wave propagation in heterogeneous elastic and anisotropic models with TTI or TORT symmetry, for example. This new method can be naturally extended to deal, for example, with tilted Monoclinic symmetry cases as well, by simply using elastic tensors with monoclinic symmetry, which are characterized by 13 independent non-zero elements, while these are five for TTI and nine for TORT. By conceptualizing velocity models as frame bundles, where local coordinate frames align with the natural frame of elastic tensors, this new method generalizes and improves upon current practice of using directional derivatives. This new method encodes frame rotations directly into covariant wave equations, ensuring differential operators are compatible with the metric of the new manifold formed by combining spatial dimensions with the orientation manifold. This new method may enhance stability and / or reduce numerical artifacts in wavefield simulations, potentially improving reflectivity and velocity model estimation in inverse problems like LSRTM or FWI. This new method also provides several insights and benefits for wave propagation in anisotropic media, for example: it separates the effects of the underlying space geometry from those of elastic tensor orientation, lends physical meaning to orientation changes of symmetry elements as new tidal forces, introduces a new model parameter constraint in the form of a curvature tensor for inversion problems, and opens avenues for developing more accurate and efficient solvers that can adaptively account for frame rotations. Additionally, this new method encourages further investigation into orientation parameterization, such as using Euler angles or quaternions.

[0130] Figure 3 schematically depicts a method 300 according to an implementation for simulating seismic wave propagation in media. The method is implemented by a computer comprising a processor and a memory. At step 302, the method comprises obtaining an earth model having a grid. At step 304, the method 300 includes modelling, using the earth model, a medium as a set of elastic tensors, wherein respective elastic tensors of the set thereof have tilted orientations relative to the grid. At step 306, the method 300 simulating seismic wave propagation in the medium using a waveequation operator, wherein the wave-equation operator is invariant to orientation changes of the respective elastic tensors of the set thereof. At step 308, the method 300 includes outputting a result of simulating seismic wave propagation, such as for example using a display.

[0131] It is to be understood that the specific order or hierarchy of operations in the method depicted in FIG. 3 and throughout this disclosure are instances of example approaches and can be rearranged while remaining within the disclosed subject matter. For instance, any of the operations depicted in FIG. 3 may be omitted, repeated, performed in parallel, performed in a different order, and / or combined with any other of the operations depicted in FIG. 3 or discussed herein.

[0132] Turning to FIG. 4, a system 400 to process communication data can include one or more computing devices 402 for performing the techniques discussed herein. In one implementation, the one or more computing devices 402 include a computing device and / or one or more servers of a processing system to generate and execute the techniques discussed herein as a software application and / or a module or algorithmic component of software.

[0133] In some instances, the computing device 402 can include a computer, a personal computer, a desktop computer, a laptop computer, a terminal, a workstation, a server device, a cellular or mobile phone, a mobile device, a smart mobile device a tablet, a wearable device (e.g., a smart watch, smart glasses, a smart epidermal device, etc.) a multimedia console, a television, an Internet-of-Things (loT) device, a smart home device, a medical device, a virtual reality (VR) or augmented reality (AR) device, a vehicle (e.g., a smart bicycle, an automobile computer, etc.), and / or the like. The computing device 402 may be integrated with, form a part of, or otherwise be associated with the system 400. It will be appreciated that specific implementations of these devices may be of differing possible specific computing architectures not all of which are specifically discussed herein but will be understood by those of ordinary skill in the art.

[0134] The computing device 402 may be a computing system capable of executing a computer program product to execute a computer process. Data and program files may be input to the computing device 402, which reads the files and executes the programs therein. Some of theelements of the computing device 402 include one or more processors 404, one or more memory devices 406, and / or one or more ports, such as input / output (IO) port(s) 408 and communication port(s) 410. Additionally, other elements that will be recognized by those skilled in the art may be included in the computing device 402 but are not explicitly depicted in FIG. 4 or discussed further herein. Various elements of the computing device 402 may communicate with one another by way of the communication port(s) 410 and / or one or more communication buses, point-to-point communication paths, or other communication means.

[0135] The processor 404 may include, for example, a central processing unit (CPU), a microprocessor, a microcontroller, a digital signal processor (DSP), and / or one or more internal levels of cache. There may be one or more processors 404, such that the processor 404 comprises a single central-processing unit, or a plurality of processing units capable of executing instructions and performing operations in parallel with each other, commonly referred to as a parallel processing environment.

[0136] The computing device 402 may be a conventional computer, a distributed computer, or any other type of computer, such as one or more external computers made available via a cloud computing architecture. The presently described technology is optionally implemented in software stored on the data storage device(s) such as the memory device(s) 406, and / or communicated via one or more of the I / O port(s) 408 and the communication port(s) 410, thereby transforming the computing device 402 in FIG. 4 to a special purpose machine for implementing the operations described herein. Moreover, the computing device 402, as implemented in the system 400, receives various types of input data (e.g., the input data, well data, etc.) and transforms the input data through various stages of the data flow into new types of data files (e.g., the simulation). Moreover, these new data files are transformed further into the prediction data and sent to the computing device 104 to provide information regarding the simulation, which enables the computing device 402 to do something it could not do before — simulating seismic wave propagation in media.

[0137] Additionally, the systems and operations disclosed herein represent an improvement to the technical field of machine modeling. For instance, the system 400 can generate estimate data with less data and fewer computing resources without human intervention. Moreover, data can be leveraged to provide a highly efficient and effective estimation of data. These techniques are rooted in technology and could not have existed prior to the advent of simulation software.

[0138] The one or more memory device(s) 406 may include any non-volatile data storage device capable of storing data generated or employed within the computing device 402, such as computerexecutable instructions for performing a computer process, which may include instructions of both application programs and an operating system (OS) that manages the various components of the computing device 402. The memory device(s) 406 may include, without limitation, magnetic disk drives, optical disk drives, solid state drives (SSDs), flash drives, and the like. The memory device(s) 406 may include removable data storage media, non-removable data storage media, and / or external storage devices made available via a wired or wireless network architecture with such computer program products, including one or more database management products, web server products, application server products, and / or other additional software components. Examples of removable data storage media include Compact Disc Read-Only Memory (CD-ROM), Digital Versatile Disc Read-Only Memory (DVD-ROM), magneto-optical disks, flash drives, and the like. Examples of nonremovable data storage media include internal magnetic hard disks, SSDs, and the like. The one or more memory device(s) 406 may include volatile memory (e.g., dynamic random access memory (DRAM), static random access memory (SRAM), etc.) and / or non-volatile memory (e.g., read-only memory (ROM), flash memory, etc.).

[0139] Computer program products containing mechanisms to effectuate the systems and methods in accordance with the presently described technology may reside in the memory device(s) 406 which may be referred to as machine-readable media. It will be appreciated that machine-readable media may include any tangible non-transitory medium that is capable of storing or encoding instructions to perform any one or more of the operations of the present disclosure for execution by a machine or that is capable of storing or encoding data structures and / or modules utilized by or associated with such instructions. Machine-readable media may include a single medium or multiple media (e.g., a centralized or distributed database, and / or associated caches and servers) that store the one or more executable instructions or data structures.

[0140] In some implementations, the computing device 402 includes one or more ports, such as the I / O port(s) 408 and the communication port(s) 410, for communicating with other computing, network, or vehicle computing devices. It will be appreciated that the I / O port 408 and the communication port 410 may be combined or separate and that more or fewer ports may be included in the computing device 402.

[0141] The I / O port 408 may be connected to an I / O device, or other device, by which information is input to or output from the computing device 402. Such I / O devices may include, without limitation, one or more input devices, output devices, and / or environment transducer devices.

[0142] In one implementation, the input devices convert a human-generated signal, such as, human voice, physical movement, physical touch or pressure, and / or the like, into electrical signals as input data into the computing device 402 via the I / O port 408. Similarly, the output devices may convert electrical signals received from the computing device 402 via the I / O port 408 into signals that may be sensed as output by a human, such as sound, light, and / or touch. The input device may be an alphanumeric input device, including alphanumeric and other keys for communicating information and / or command selections to the processor 404 via the I / O port 408. The input device may be another type of user input device including, but not limited to direction and selection control devices, such as a mouse, a trackball, cursor direction keys, a joystick, and / or a wheel; one or more sensors, such as a camera, a microphone, a positional sensor, an orientation sensor, an inertial sensor, and / or an accelerometer; and / or a touch-sensitive display screen (“touchscreen”). The output devices may include, without limitation, a display, a touchscreen, a speaker, a tactile and / or haptic output device, and / or the like. In some implementations, the input device and the output device may be the same device, for example, in the case of a touchscreen.

[0143] The environment transducer devices convert one form of energy or signal into another for input into or output from the computing device 402 via the I / O port 408. For example, an electrical signal generated within the computing device 402 may be converted to another type of signal, and / or vice-versa. In one implementation, the environment transducer devices sense characteristics or aspects of an environment local to or remote from the computing device 402, such as, light, sound, temperature, pressure, magnetic field, electric field, chemical properties, physical movement, orientation, acceleration, gravity, and / or the like.

[0144] In one implementation, the communication port 410 is connected to the network(s) 112 so the computing device 402 can receive network data useful in executing the methods and systems set out herein as well as transmitting information and network configuration changes determined thereby. Stated differently, the communication port 410 connects the computing device 402 to one or more communication interface devices configured to transmit and / or receive information between the computing device 402 and other devices by way of one or more wired or wireless communication networks or connections. Examples of such networks or connections include, without limitation, Universal Serial Bus (USB), Ethernet, Wi-Fi, Bluetooth®, Near Field Communication (NFC), and so on. One or more such communication interface devices may be utilized via the communication port 410 to communicate with one or more other machines, either directly over a point-to-point communication path, over a wide area network (WAN) (e.g., the Internet), over a local area network (LAN), over a cellular network (e.g., third generation (3G), fourth generation (4G), Long-Term Evolution (LTE), fifth generation (5G), etc.) or over another communication means. Further, thecommunication port 410 may communicate with an antenna or other link for electromagnetic signal transmission and / or reception.

[0145] In an example, software, modules, services, and operations discussed herein may be embodied by instructions stored on the memory device(s) 406 and executed by the processor 404.

[0146] The system set forth in FIG. 4 is but one possible example of a computing device 402 or computer system that may be configured in accordance with aspects of the present disclosure. It will be appreciated that other non-transitory tangible computer-readable storage media storing computer-executable instructions for implementing the presently disclosed technology on a computing system may be utilized. In the present disclosure, the methods disclosed may be implemented as sets of instructions or software readable by the computing device 402.

[0147] Although exemplary implementations have been shown and described, it will be appreciated by those skilled in the art that various changes and modifications might be made without departing from the scope of the disclosure, as defined in the appended claims and as described above.

[0148] At least some of the implementations described herein may be constructed, partially or wholly, using dedicated special-purpose hardware. Terms such as ‘component’, ‘module’ or ‘unit’ used herein may include, but are not limited to, a hardware device, such as circuitry in the form of discrete or integrated components, a Field Programmable Gate Array (FPGA) or Application Specific Integrated Circuit (ASIC), which performs certain tasks or provides the associated functionality. In some embodiments, the described elements may be configured to reside on a tangible, persistent, addressable storage medium and may be configured to execute on one or more processors. These functional elements may in some embodiments include, by way of example, components, such as software components, object-oriented software components, class components and task components, processes, functions, attributes, procedures, subroutines, segments of program code, drivers, firmware, microcode, circuitry, data, databases, data structures, tables, arrays, and variables. Although the example embodiments have been described with reference to the components, modules and units discussed herein, such functional elements may be combined into fewer elements or separated into additional elements. Various combinations of optional features have been described herein, and it will be appreciated that described features may be combined in any suitable combination. In particular, the features of any one example embodiment may be combined with features of any other embodiment, as appropriate, except where such combinations are mutually exclusive. Throughout this specification, the term “comprising” or “comprises” means including the component(s) specified but not to the exclusion of the presence of others.Working example - proof of concept.

[0149] We will in the following present a numerical experiment on a simple 2D velocity model, demonstrating that the orientation-invariant wave equations described herein can match or outperform traditional methods.

[0150] In the example, the model consists of a single homogeneous block of transversely isotropic (Tl) material. All parameters are fixed except for the tilt of the Tl symmetry axis, which jumps from 0 to 90 degrees midway through the model (at a depth of around 500 m): the upper half has a vertical symmetry axis (VTI), and the lower half a horizontal symmetry axis (HTI), as shown in Figure 5. This orientation jump creates an interface, introducing an impedance contrast that generates reflections and mode conversions.For reference, Table 1 summarizes the model and simulation parameters.Table 1 - Model and simulation parameters

[0151] Parameter

[0152] Symbol

[0153] Value

[0154] Density

[0155] p

[0156] 1.0 g / cm3

[0157] P-wave velocity

[0158] Vp

[0159] 3.6 km / s

[0160] S-wave velocity

[0161] K;

[0162] 1.8 km / s

[0163] Thomsen parameter

[0164] 6

[0165] 0.23

[0166] Thomsen parameter

[0167] 6

[0168] 0.17

[0169] Grid spacing

[0170] Ax = Az

[0171] 10 m

[0172] Source location

[0173] (x,z)

[0174] (1,000 m, 10 m

[0175] Source wavelet

[0176] -

[0177] Ricker, 10 Hz

[0178] We propagate waves using the following second-order wave equation for the stress tensor < Jij (From Equation 3):dt^ij = ~A^A0kl,

[0179] simulated numerically with finite differences in time domain. For the example, two approaches are compared for evaluating the Laplacian A:1. Pre-rotation (traditional)

[0180] Rotate the density-normalized elastic tensor A1to the global frame, then compute the Euclidean Laplacian with three Cartesian coordinates:y1d2^ ■= £ / —i ~ drx~fr.k=lK2. Orientation-invariant (novel approach)

[0181] Keep A1- - in the local frame and use the presented approach to account for orientation changes directly by including the orientation of the elastic tensors as addional coordinates.

[0182] Figure 6 shows how the orientation jump at 500 m depth affects the elastic stiffness components in the global frame (traditional approach), establishing a mental picture of the subsurface model in which the stress wavefields will propagate. Henceforth, we refer to the density-normalized elastic coefficients A - using the more convenient Voigt notation, which condenses indices pairs into a single index, like so: 11 -> 1; 13 -> 5; 33 -» 3. To make the Voigt notation change more explicit we also modify the letter denoting the tensor from A to C. Thus A^ = C1X, A™ =Ai3 =cis> andso onandsoforth.

[0183] In Figure 6, C33and Ctlswap values across the interface, creating the expected impedance contrast. The remaining components contributing to wave propagation in the example (C13, C15> f-35 > ^5s) should be invariant to the 90-degree rotation; nevertheless, the smoothing of the orientation jump, done to stabilize the finite difference scheme, has the unintended side effect of adding a contrast in those components as well. This is a common pitfall of finite differences in general and the first approach in particular: strong contrasts often require smoothing, but this may introduce false reflectors / scatterers which in turn produce mode conversion events that should not exist. The orientation-invariant approach is less affected by such smoothing, because it keeps the elastic tensor in the local frame and handles orientation changes through the covariant differential operators.

[0184] Figures 7 and 8 show snapshots at t = 225 ms of the three stress components from simulations done for both approaches.

[0185] Figure 7 shows snapshots of stress wavefield generated by a pressure source for horizontal, shear, and vertical components at t=225 ms (pre-rotation approach). Symmetry axis changes induce PF, PS, and SP mode conversions.

[0186] Figure 8 shows snapshots of stress wavefields at t=225 ms from simulations with covariant Laplacian (orientation-invariant approach), which also produces PP-wave reflection and PS and SP mode conversions from the change in symmetry axis orientation.

[0187] An explosive source generates P- and S-wavefronts, which interact with the orientation interface at 500 m, producing PP, PS, and SP reflections and transmissions as expected. To wit, P-wavefronts have the same sign across the normal stress components au; S-wavefronts have opposite signs in those normal stress components or appear in the shear stress one. Both approaches yield similar kinematic behavior, and the bulk of the differences come from amplitude mismatches caused by subtly distinct energy partitioning at the model interfaces and boundaries.

[0188] Note that in the orientation-invariant approach, wavefields are in the local frame, while in the pre-rotation approach, they are in the global frame. Thus, comparing normal stress responses in the lower half of the model requires swapping horizontal and vertical components between Figures 7 and 8. Essentially, orientation-invariant approach reproduces the results from pre-rotation approach, but without rotating the elastic tensor to the global frame.

[0189] Figure 9 shows the snapshots of the shear-stress component wavefields at t=225 ms from simulations with pre-rotation (traditional) approach (left) and proposed orientation invariant approach (right), but with boosted amplitude.

[0190] Differences in the shear stress can be examined directly (because of the invariance of a13to the 90-degree rotation) and they showcase how each approach handles the orientation jump. In the pre-rotation approach, incorrect contrasts in off-diagonal stiffness and in C55cause less coherent and smaller transmitted wavefields from the incident shear-wave energy across the interface because of SS- and SP-wave reflections.

[0191] The orientation-invariant approach, on the other hand, shows much stronger coupling between normal and shear stress components, with a more coherent shear wavefront crossing the interface without conspicuous reflections. The phases of reflected and transmitted PS-wavefields are opposite to each other in the orientation-invariant approach, as expected, which is not the case in the pre-rotation approach. Furthermore, at the bottom of the model, the incident P-wave creates reflection responses only in normal components in the pre-rotation approach, while in the orientation-invariant approach we observe responses in all three components, again as expected.

[0192] To remove the coordinate frame effect from the comparison, we show principal stress wavefields (i.e., the eigenvalues of the 2D stress tensor) for both methods in Figures 10 (prerotation) and 11 (orientation-invariant), along with their differences (Figure 12). The principal stress field snapshots prove that the two methods yield comparable results: the norm of the difference is less than 3% of the norm of each individual principal stress wavefield and differences becomevisible only when amplitudes are boosted tenfold or more, as in Figure 12. Furthermore, this figure confirms that the bulk of the differences stem from energy partitioning at the interface and at the bottom edge, documenting that the pre-rotation approach tends to produce stronger PP-events at the expense of SS- and PS- reflections and transmissions. The orientation-invariant approach yields more consistent amplitudes and sign changes, better representing the true physical behavior.

[0193] In summary, the presented orientation-invariant method not only matches traditional results but surpasses them: covariant differential operators are more robust to orientation changes, avoiding errors from rotation and smoothing of elastic tensor components required for stability in the traditional approach.

[0194] In the present disclosure, any term of degree such as, but not limited to, “substantially,” as used in the description and the appended claims, should be understood to include an exact, or a similar, but not exact configuration. Similarly, the terms “about” or “approximately,” as used in the description and the appended claims, should be understood to include the recited values or a value that is three times greater or one third of the recited values. For example, about 3 mm includes all values from 1 mm to 9 mm, and approximately 50 degrees includes all values from 16.6 degrees to 150 degrees.

[0195] Lastly, the terms “or” and “and / or,” as used herein, are to be interpreted as inclusive or meaning any one or any combination. Therefore, “A, B, or C” or “A, B, and / or C” mean any of the following: “A,” “B,” or“C”; “A and B”; “A and C”; “B and C”; “A, B and C.” An exception to this definition will occur only when a combination of elements, functions, steps or acts are in some way inherently mutually exclusive.

[0196] While the present disclosure has been described with reference to various implementations, it will be understood that these implementations are illustrative and that the scope of the present disclosure is not limited to them. Many variations, modifications, additions, and improvements are possible. More generally, implementations in accordance with the present disclosure have been described in the context of particular implementations. Functionality may be separated or combined differently in various implementations of the disclosure or described with different terminology. These and other variations, modifications, additions, and improvements may fall within the scope of the disclosure as defined in the claims that follow.

Claims

CLAIMS1. A method of simulating seismic wave propagation in media, the method implemented by a computer comprising a processor and a memory, the method comprising:modelling, using an earth model having a grid, a medium as a set of elastic tensors, wherein respective elastic tensors of the set thereof have tilted orientations relative to the grid; and wherein tensor orientation is added to the coordinates of the grid so that the model comprises three cartesian coordinates and three local elastic tensor orientation coordinates; andsimulating seismic wave propagation in the medium using a wave-equation operator, wherein the wave-equation operator is invariant to orientation changes of the respective elastic tensors of the set thereof.

2. The method according to claim 1, wherein simulating seismic wave propagation in the medium using the wave-equation operator comprises aligning respective local frames of the wave-equation operator with respective symmetry elements of the respective elastic tensors of the set thereof.

3. The method according to claim 2, wherein aligning respective local frames of the wave-equation operator with the respective symmetry elements comprises providing the respective local frames by a fiber bundle.

4. The method according to claim 3, wherein providing the respective local frames by the fiber bundle comprises attaching a copy of an orientation grid (for example, of the respective symmetry elements) to respective points in a spatial coordinate grid of the earth model.

5. The method according to claim 3, wherein providing the respective local frames by the fiber bundle comprises representing space of all orientations of respective local frames using the respective fibers.

6. The method according to claim 3, wherein providing the respective local frames by the fiber bundle comprises representing the respective local frames axially and complementary to the respective fibers.

7. The method according to claim 3, wherein providing the respective local frames by the fiber bundle comprises providing the respective local frames by the fiber bundle having six dimensions, for example generally curved manifolds; wherein a base space provides coordinates for location in 3Dspace; and wherein respective fibers represent orientation manifolds, for example parameterized by coordinates and describing respective local orientations of the respective symmetry elements.

8. The method according to one of the previous claims, wherein the medium has Tilted Transverse Isotropy, TTI, symmetry, Tilted Orthorhombic, TORT, symmetry and / or Tilted Monoclinic symmetry.

9. The method according to claim 1, wherein the earth model comprises and / or is a heterogeneous and / or anisotropic earth model.

10. The method according to claim 1, comprising outputting a result of simulating seismic wave propagation in the medium.

11. A method of seismic imaging of an earth model, comprising the method according to claim 1.

12. A method of seismic inference of an earth model, comprising the method according to claim 1.

13. A method of adjusting an earth model, comprising the method according to claim 1.

14. A method of inversion of seismic data, comprising the method according to claim 1.

15. A computer, comprising a processor and a memory, configured to perform the method according to any of claims 1 to 14;a computer program comprising instructions which, when executed by a computer, comprising a processor and a memory, cause the computer to perform the method according to any of claims 1 to 14; and / ora tangible (optionally, a non-transitory) computer-readable recording medium storing instructions which, when executed by a computer, comprising a processor and a memory, cause the computer to perform the method according to any of claims 1 to 14.

Images