A Joint Inversion-Segmentation approach to Assisted Seismic Interpretation
AA J
OINT I NVERSION -S EGMENTATION APPROACH TO A SSISTED S EISMIC I NTERPRETATION
A P
REPRINT
Matteo Ravasi
KAUSTThuwal, Kingdom of Saudi Arabia [email protected]
Claire Birnie
KAUSTThuwal, Kingdom of Saudi Arabia [email protected]
February 9, 2021 A BSTRACT
Structural seismic interpretation and quantitative characterization are historically intertwined pro-cesses. The latter provides estimates of properties of the subsurface which can be used to aidstructural interpretation alongside the original seismic data and a number of other seismic attributes.In this work, we redefine this process as a inverse problem which tries to jointly estimate subsurfaceproperties (i.e., acoustic impedance) and a piece-wise segmented representation of the subsurfacebased on user-defined macro-classes . By inverting for the quantities simultaneously, the inversionis primed with prior knowledge about the regions of interest, whilst at the same time it constrainsthis belief with the actual seismic measurements. As the proposed functional is separable in the twoquantities, these are optimized in an alternating fashion, where each subproblem is solved usinga Primal-Dual algorithm. Subsequently, each class is used as input to a workflow which aims toextract the perimeter of the detected shapes and to produce unique horizons. The effectiveness of theproposed method is illustrated through numerical examples on synthetic and field datasets.
Seismic interpretation forms part of an integrated workflow for reservoir mapping and characterization, and representsone of the most important steps for the successful exploration of underground resources and management of existingones. It is traditionally divided into two components: structural and quantitative interpretation. The former aims atidentifying the geological framework (i.e., horizons and faults) whilst the latter focuses on retrieving properties of thesubsurface.Whilst structural interpretation is nowadays mostly carried out manually and led by geological understanding, severalkey geophysical principles are also followed early on to identify an interpretation strategy to guide to the interpreter’sjob: with the aid of sonic and density logs, a first step in an interpretation project is represented by the creation ofwell-ties, which allow interpreters to link the geological stratigraphy to the observed seismic response - in other words,interpreters look for changes in acoustic and shear impedances (or velocity ratios) in log data, which trigger a seismicresponse, and perform seismic modelling to identify such events in the observed data. However, as the word implies, interpretation is a subjective matter whose ultimate goal is that of reconstructing the geological story contained in aseismic volume [1].With recent advances in artificial intelligence and deep learning, several attempts at automatic seismic interpretationhave been reported [2, 3, 4, 5]. A common feature of most of these approaches is that of turning a geologically-drivenprocess into a computer vision and pattern recognition problem. They have been shown to be very successful for somespecific tasks where a clear pattern can be observed and learned - e.g., for top salt interpretation [2], where a seismicrich area is abruptly interrupted by a strong reflector, followed by an area with little to no seismic signal). On the otherhand, these methods currently fail at data integration, something that human interpreters are very good at. The use ofwell-ties is a good example of this - where well-tie principles are not directly encoded in the machine learning trainingprocess resulting in the need for interpreters inputs prior to training. Typically, interpreters are still required to define the a r X i v : . [ phy s i c s . g e o - ph ] F e b oint Inversion-Segmentation for Seismic Interpretation A PREPRINT horizons of interest and in some cases interpret it for a number of lines in order to train a machine that is subsequentlyable to auto-track the same horizons of interest in the rest of the dataset [3].Post- and pre-stack inversion represent the workhorses in the process of quantitative characterization of the subsurface.Whilst coming in many different flavors, most commonly used methods estimate parameters by optimizing a least-squares functional with spatial regularization terms and, in the pre-stack case, extra regularization terms aimed atenforcing correlation between parameters [6]. However, since the input data is band-limited in nature and lacks both thevery low and very high frequencies, least-squares approaches suffer from the inability of the inversion to fully recoverthe entire spectrum, leading to unwanted oscillatory behaviors near interfaces in the recovered model. While a properlydefined low frequency background model can partially fill such a frequency gap, the remaining missing frequenciescan only be compensated by providing additional prior information in the form of more suitable regularizers. TVregularization is one such type of regularization as it favours blockiness in the recovered impedance model; addingsuch regularization does however come at the cost of making the overall functional non-smooth and when used alone isknown to lead to a systematic loss of contrast [7]. The Split-Bregman algorithm [8] has been successfully employed insolving L1-regularized problems in a variety of fields and it is used by [9] in the context of seismic impedance inversionof reflectivity data (after a preliminary step of blind deconvolution). Similarly, [10] and [11] use the AlternatingDirection Method of Multipliers (ADMM) for both post- and pre-stack inversion of band-limited data.In this work, we propose a different approach to computer-aided seismic interpretation, which aims to reconcile theclassical principles of seismic interpretation to the underlying optimization problem that we define and solve. Theinterpreter input is here limited to the definition of acoustic zones identified by a single acoustic impedance (AI) value(or elastic zones identified by triplets of parameters) that we wish to ultimately extract from the seismic data. A jointinversion and segmentation of the input seismic data is then performed and both inverted acoustic (or elastic) propertiesand seismic horizons are produced. Examples on synthetic and field datasets show that our methodology can successfullyinterpret key structures as defined by the interpreter in a fully automatic manner without any manual labelling orexpensive upfront training. Moreover, as a by-product, both an acoustic impedance model and a segmentation modelare retrieved that can be used to further refine the interpreted horizons as well as inputs to subsequent steps of reservoirmodelling.
Traditionally, seismic data has been interpreted by looking at a variety of attributes (e.g., full stack, partial stacks,seismic attributes) and most commercial softwares are built around this principle. Whilst used in some cases to aidinterpretation, elastic properties derived from seismic pre- or post-stack inversion are generally not considered themain source of information against which to interpret horizons. Nevertheless, these properties are the physical reasonbehind what is observed within the seismic data: a large increase (or decrease) in acoustic impedance at a geologicalboundary produces a strong seismic response; conversely, a small change cannot trigger a response and so our ability tointerpret such geological boundary is in vain. We therefore ask ourselves whether a machine, unhindered by humanhabits or limitations of current softwares, could be more successful at aiding interpretation if working directly withelastic properties other than the original seismic data? Or in other words, can we turn the problem of interpreting anhorizon from following wiggles on seismic sections into a segmentation task for the acoustic (or elastic) properties inthe geological formations of interest? In the next section we propose a mathematical formulation that translates intoequations our objective.
Mathematical formulation
Inspired by the work of [12] on joint reconstruction and segmentation in the context of medical imaging, we definea functional to jointly invert seismic data for their acoustic impedance as well as to produce a segmentation of thesubsurface into N c classes, each of which is characterized by a different range of impedance values. The underlyingidea of inverting for these two parameters simultaneously is to inform the inversion with prior knowledge about theregions of interest, whilst at the same time constraining this belief with the seismic measurements. The functional forour specific problem can be written as: ( m , V ) = argmin m , V ∈ C || d − Gm || + αT V ( m ) + δ N c (cid:88) j =1 N x N z (cid:88) i =1 V ji ( m i − c j ) + β N c (cid:88) j =1 T V ( V Tj ) (1)where m is a vector of size N x N z × that contains the natural logarithm of the acoustic impedance values in the areaof the subsurface of interest, V is a matrix of size N c × N x N z whose columns contain the probability of each point inthe subsurface to belong to each class, and c is a vector of size N c × which contains the acoustic impedance values2oint Inversion-Segmentation for Seismic Interpretation A PREPRINT associated to each class. G = WD is the post-stack modelling operator composed of a first derivative operator D and aconvolution operator W whose convolution kernel is the estimated wavelet w divided by 2, and d is post-stack seismicdata of size size N x N z × . Here, we use the convention that x i is the i-th element of a vector and X i represents theextraction of the i-th column of a matrix (whilst X Tj is the j-th row of a matrix transposed into a column vector).Furthermore, we further constrain each column of V to be in the unit Simplex, C = { V i ∈ R + : (cid:80) N c j =1 V ji = 1 } ∀ i =1 , , ..., N z N x . Note that by ensuring that the sum of the elements of each column of V is equal to 1 and every elementis positive, we can interpret such values as the probability of each point in the subsurface belonging to a certainclass. The Total Variation (TV) regularization term is defined as T V ( x ) = ||∇ x || , = (cid:80) N x N z i =1 (cid:112) ( x x ) i + ( x z ) i where ∇ : R N x N z → R × N x N z is the gradient operator that transforms a vector into a matrix whose rows contain thederivatives ( x x ) and ( x z ) computed along the x- and z directions, respectively. Finally, α , β , and δ are regularizationparameters used to balance the inversion and segmentation terms of the functional. In order to solve the functional in equation 1 we note that, whilst non-convex in the joint argument ( m , V ) , this costfunction is convex in each individual variable space. A splitting approach is therefore used to minimise the two convexproblems in an alternating fashion. Moreover, since solving TV-regularized problems can lead to a systematic loss ofcontrast in the estimated model [13], we replace both TV regularization terms with their generalized Bregman distance( D p T V ( u , u (cid:48) ) = T V ( u ) − T V ( u (cid:48) ) − ( u − u (cid:48) ) T p ) and Bregman iterations are introduced to solve each independentminimisation problem [7]. The overall algorithm reads as: m k = argmin m || d − Gm || + α ( T V ( m ) − m T p k − ) + δ N c (cid:88) j =1 N x N z (cid:88) i =1 V k − ji ( m i − c j ) (2a) p k = p k − − α G T ( Gm k − d ) + 2 δ N c (cid:88) j =1 V Tj ( m k − c j ) (2b) V k = argmin V ∈ C δ N c (cid:88) j =1 N x N z (cid:88) i =1 V ji ( m ki − c j ) + β N c (cid:88) j =1 ( T V ( V Tj ) − ( V Tj ) T Q T k − j ) (2c) Q kji = Q k − ji − δβ ( m ki − c j ) ∀ j = 1 , , .., N c i = 1 , , .., N x N z (2d)where p and Q are the sub-gradients of their corresponding Bregman distances. In the following, we will focus on thesolution of equations 2a and 2c. The Primal-Dual solver
Proximal algorithms are a family of solvers able to minimize non-smooth, convex functionals like those in equations 2aand 2c [14]. Minimisation is performed in an iterative fashion and each iteration relies on the ability of evaluating theso-called proximal operator of the functional f to be minimised: prox τf ( u ) = argmin x f ( x ) + 12 τ || x − u || (3)As we will see in the following, closed form solutions exist for many commonly used functionals, meaning that thesesub-problems can be solved very quickly with fast specialized methods.The Primal-Dual algorithm [15] is a special type of proximal solver that shows O (1 /N ) rate of convergence in finitedimensions. It works by recasting any functional of this kind: argmin x f ( Kx ) + (cid:88) i x T z i + g ( x ) (4)into its primal-dual equivalent (a saddle-point problem): argmin x argmax y y T ( Kx ) + (cid:88) i z Ti x + g ( x ) − f ∗ ( y ) (5)where f and g are convex (possibly non-smooth) functionals, f ∗ is the convex conjugate of f and K is a linear operatorthat maps the vector x ∈ R m into a vector y ∈ R n . 3oint Inversion-Segmentation for Seismic Interpretation A PREPRINT
A series of iterations is then introduced in order to obtain convergence at the saddle point: y k +1 = prox µg ∗ ( y k + µ K ¯ x k ) x k +1 = prox τf ( x k − τ ( K H y k +1 + (cid:80) i z i ))¯ x k +1 = x k +1 + θ ( x k +1 − x k ) (6)where θ ∈ [0 , ( θ = 1 will be used in our work), τ and µ represent the step-lengths of the two sub-gradients. Toensure convergence τ µL < ( L = || K || = λ max ( K H K ) - i.e., the spectral radius of the linear operator. We start by analysing equation 2a and showing that it is possible to write it to be on the form of equation 4. Specifically,we write K = ∇ , z = − α p k − , f = || || , , and g = || d − Gm || + δ (cid:80) N c j =1 (cid:80) N x N z i =1 V k − ji ( m i − c j ) . Forsimplicity, we rearrange some of the terms in g such that we can write it as g = || ˆ d − ˆ Gm || : ˆ V = diag { (cid:112) V T } ...diag { (cid:113) V TN c } , ˆ c = (cid:112) V T c ... (cid:113) V TN c c N c , ˆ G = (cid:20) G √ δ ˆ V (cid:21) , ˆ d = (cid:20) d √ δ ˆ c (cid:21) (7)To conclude the proximal operator of g is simply that of a quadratic functional which can be written as: prox τg ( m ) = (cid:16) I + τ ˆ G T ˆ G (cid:17) − (cid:16) m + τ ˆ G T ˆ d (cid:17) (8)Similarly, the proximal operator of f is that of a L , norm - i.e. sum of the Euclidean norms of the columns of a matrix Y - which is: prox τf ( Y ) = ( prox τf (( Y T ) ) T , ..., prox τf (( Y T ) N ) T ) T , prox τf ( y ) = (cid:18) − στmax {|| y || , τ } (cid:19) y (9) To start we rewrite equation 2c as follows: v k = argmin v ∈ C v T ( δ g − β q k − ) + β N c (cid:88) j =1 T V ( V Tj ) (10)where v = V ec ( V ) and q = V ec ( Q ) are the vectorized version of V and Q , respectively, obtained by concatenatingtheir rows into a vector. Moreover, we define the weighting vector g = (( m − c ) , .., ( m − c N c ) ) T , where thesuperscript is used to indicate the squared difference between each element of the vector m and the chosen coefficient c j . In its current form, equation 10 can be shown to be on the form of equation 4. This is achieved by writing K = ∇ , z = δ g − β q k − , f = || || , , and g = i C i.e., the indicator of the Simplex function.The algorithm described in equation 2 is implemented here using the PyLops framework [16]. Horizon extraction
Horizon extraction is performed per class as illustrated in Figure 1. First, a binary class image is computed from thesegmented image identifying the areas which are allocated to that class after the segmentation is complete. The binaryclass image undergoes two cleaning steps: a feature selection procedure that removes small closed foreground objectsand morphological erosion which converts foreground pixels to the background pixel value when the background isdominant within a fixed window size. The terminology of foreground and background here represents the areas in thebinary class image that have either been detected as that class or not.Subsequently, the previously defined Total Variation is computed on the cleaned class image. Such a measure isfrequently used for edge detection algorithms as it enhances the edges of objects isotropically. The TV computationidentifies the boundaries of areas belonging to that class, i.e., the edges of the foreground objects, such that pixels witha high TV value are indicative of a boundary. Finally, a threshold is applied where every pixel above the threshold isdetermined as a potential horizon point. 4oint Inversion-Segmentation for Seismic Interpretation A PREPRINT
Once detected, neighbouring horizon points are joined in a left-right, top-bottom fashion. Beginning at the top,most-left horizon point, ( x , y ) , neighbouring horizon points are identified within the limits x < x < x + 1 and y − < y < y + 1 . If there are multiple neighbouring points satisfying the TV threshold, the one with the highestTV value is selected as the next horizon point and the neighbour search begins again with the new points, ( x , y ) ,searching the space, x < x < x + 1 and y − < y < y + 1 . The neighbour search continues until no points aredetected in the neighbour space, concluding that line. At this point, a new neighbour search begins from the highest,most-left remaining horizon points to create another line. The colours and separated lines in the connect image of Figure1 represent the results of the neighbour search step. Once the neighbour search is exhausted, the lines are labelled usingthe segmented image to determine the class above and class below the lines, as illustrated in the label image of Figure1. Three options exist for labelling the lines depending on the complexity of the geology and the homogeneity of themedium above and below the line. The labelling options are: only consider the class above, only consider the classbelow, or consider both classes above and below.Similar to the connect step performed on a pixel level, a secondary combining procedure is performed on the labelledlines considering any lines that begin within a specified window from the ending of the previous line. Due to thesmaller number of lines available to be joined, in comparison to the pixels in the connect ’ step, a larger window can beconsidered allowing joining of points previously too far apart, as illustrated by the combining of the 3 AB lines. Forlines that are too far apart, such as BC , these will remain as separate lines to avoid any artefacts being generated in theregridding step.Finally, the regridding step is performed for each identified horizon line. For continuous horizons, such as AB , theregridding procedure ensures a single depth value per horizontal location. On the other hand, when horizons presentgaps (e.g., BC ), the regridded horizon is further masked to avoid the creation of spurios connections. As a result, all ofthe identified horizons span the full lateral extent of the model.Once the horizons have been labelled and extracted for each class, duplicate horizons exist where they have beenidentified as the top horizon from one class image and as the bottom horizon from another class image. The TV imagevalues vary based on the class image therefore duplicate horizons are unlikely to be identical. For example in the case ofthe BC horizon, the fault may have been detected when the class image represented the medium below the BC horizon.Duplicate horizons arising from the above mentioned scenario are undesirable however there are other scenarios wherehorizons may be labelled the same based on their neighbouring classes where it is desirable to retain both, such ashighly-layered settings.A horizon selection procedure is performed as follows:1. All horizons are initially set as “to-be-kept”.2. The longest horizon, i.e. with the fewest NaN values, is identified and removed from the list of potentialhorizon lines.3. For each remaining horizon line, any shared points on the x-axis are compared and the average depth differencebetween the points is computed.4. If the difference is less than a given threshold the horizon is considered a duplicate of the longest horizon andis removed from the list of potential horizon lines.5. If the difference is larger than the threshold, the horizon line remains on the potential horizon list and is stillconsidered “to-be-kept”.6. Steps 2-5 are repeated until there are no horizons left on the potential horizon list, i.e., until all horizons linesare either selected or rejected.The selected horizons, alongside the inverted and segmented models, ultimately represent the final outputs of ouralgorithm. Examples
In this section, we apply our methodology to two synthetic examples and a 2D line of the open-source Volve dataset.
Simple model
The first example considers a very simple subsurface model composed of a number of horizontally stacked layers offsetby a normal fault. Because of its simplicity, we will show that a single step of the proposed inversion scheme suffices toobtain a satisfactory segmentation of the model into different zones and for the tracking algorithm to extract all the5oint Inversion-Segmentation for Seismic Interpretation A PREPRINT
Figure 1: Workflow for the extraction of horizon lines from a single class. The segmented image is the starting point ofthe horizon extraction procedure from which a binary class image is computed prior to cleaning and performing a TVcalculation to determine the edges of the class object. Neighbouring edge points are connected prior to being labelledbased on surrounding classes. Finally similarly labelled lines are combined before being regridded to the full horizontal,spatial sampling of the model. The blue line on the X-axis of the regridding step illustrates where horizontal values areNaN. x[m] t [ s ] (a) AI[kg/(m s)] x[m] (b) Seismic (c)
AI trace
Figure 2: a) Acoustic impedance model, b) noisy seismic data, and c) AI profile at well location (shown by the dashedblack line in panel a)horizons of interest. Nevertheless, despite its apparent simplicity, this model presents offsetted horizons which arehistorically challenging to track in seismic data.Figure 2a displays the acoustic impedance model composed of six formations with different properties. For the sake ofthis exercise, we assume that a well is drilled in the middle of the model, such that all formations are sampled. Theacoustic impedance profile is displayed in Figure 2c. The class vector c is therefore chosen to be a × vector whoseelements are the different AI values in the vertical profile. The post-stack seismic data (Figure 2b) is modelled using the G operator discussed above and a 8Hz Ricker wavelet. Note that throughout the paper we define the vertical axis interms of two-way traveltime for consistency with the field data example. Finally colored noise is added to the syntheticdataset: the noise model is created by smoothing a white noise realization along both the spatial and time axes.6oint Inversion-Segmentation for Seismic Interpretation A PREPRINT
To begin with, the noisy dataset ( d n = Gm + n ) is inverted for the underlying acoustic impedance model ( m ) using anumber of approaches:• Least-squares regularized inversion, which minimizes the following cost function: J L = || d n − Gm || + (cid:15) || m || + α ||∇ m || where ∇ is the Laplacian operator. The second term in the functional is the standardThicknov regularization, whilst the third term enforces smoothness in the solution. The model m is hereestimated by using an iterative solver for smooth, convex functionals such as LSQR using a smooth version ofthe model as starting guess. This results serves as a baseline as it represents the approach routinely used forpost-stack inversion of seismic data.• Linearized Alternating Direction Method of Multipliers (L-ADMM) with isotropic TV regularization [17]: J L − ADMM = || d n − Gm || + αT V ( m ) .• Primal-dual algorithm with isotropic TV regularization [15]: J P D = J L − ADMM
In Figure 3 the different inversion results are compared. First, we can observe that whilst computationally cheaper thanthe other methods, standard least-squares inversion suffers from a major limitation: the retrieved model shows smoothtransitions between different zones and therefore the sharp boundaries are not recovered. Whilst this is well-knownlimitation of such a method, it is particularly detrimental for our subsequent tasks of segmentation and horizon trackingas shown below. On the other hand both the L-ADMM and Primal-Dual algorithms are able to recover a model withsharp boundaries and an overall higher peak signal-to-noise ratio (
P SN R = 10 log ( N x N z max ( ˆ m ) / || m − ˆ m || )) ,where ˆ m is the estimated model). Apart from producing an estimate with slightly higher PSNR, the latter solverpossesses higher flexibility and the possibility to handle additional terms in the form of dot products as required inequation 2a of our joint inversion and segmentation algorithm that will be applied to the next example. x[m] t [ s ] (a) L2 (PSNR=52.77) x[m] (b)
L-ADMM (PSNR=54.08) x[m] (c)
PD (PSNR=55.30) (d)
AI trace
TrueBackLSL-ADMMPD
Figure 3: Inverted models via a) Least-squares regularized inversion, b) Anisotropic TV-regularized inversion Split-Bregman solver, and c) Isotropic TV-regularized inversion Primal-Dual solver. d) True, background, and inverted AIprofiles at well location.The models obtained by minimizing the first and third functionals are now used as input to the segmentation step. Inmathematical terms, this step is equivalent to solving a simplified version of equation 2c which only contains the TVnorm instead of the Bregman TV norm - in other words where Q is set to zero. Figures 4 and 5 show the results of thesegmentation for the model produced by least-squares and Primal-Dual inversions, respectively. In both cases, panel ashows the true segmented model, where each image pixel is assigned the label of its class. In panel b the binarisedversion of the estimated matrix V is displayed; binarised means that the index of the highest value of each column V is selected. The segmentation result from the least-squares model presents an evident problem inherited from thesmooth transitions across interfaces observed in the inverted model; green rings appear all around the red class due tothe fact that the model parameter of the green class lies in between that of the red and yellow classes. This becomeseven more evident when observing the probabilities of each class in the vertical pillar in the middle of the model (Figure7oint Inversion-Segmentation for Seismic Interpretation A PREPRINT x[m] t [ s ] (a) True x[m] (b)
Segmentation from LS x[m] (c)
Detected Horizons (d)
Probs.
Figure 4: a) True and b) estimated segmentation from least-squares inverted model, respectively. c) Tracked horizonsoverlaid to the noise free seismic dataset (thin lines: true horizons, thick lines: tracked horizons. d) Probabilities foreach class at well location (i.e., a single column of V at the index in the middle of the model.4d). When attempting to track horizons, the presence of such rings affect the overall quality of the estimated horizons.Moreover, some of the erroneously tracked horizons do not conform with the layering knowledge provided by the inputwell log.On the other hand, the segmentation results obtained when using the model estimated via the Primal-Dual solver aremuch more accurate. Note also that the probabilities estimated by the same number of iterations of the Primal-Dualalgorithm used for segmentation present more sharp transitions between the different zones. This greatly benefits thesubsequent step of horizon tracking which returns the five main horizons with a very high degree of precision. Overall,we can conclude that the process of identifying seismic horizons by means of tracking of the edges of closed shapedpolygons seems to be less cumbersome and more robust than the problem of tracking a horizon in a highly-oscillatorymulti-dimensional signal such as the seismic data itself. Hess model
Our second example is based on a modified version of the P-wave velocity of the SEG Hess VTI model (Figure 6a).This model is chosen for two reasons: first, we want to investigate the ability of our method to deal with stratigraphiesthat are interrupted by intrusions such as the salt body present in this model. Second, we are interested to assess theimportance of our joint inversion-segmentation approach in the presence of very strong impedance contrasts such asthose originated by the salt body.To begin with, we consider once again the case where we have complete knowledge of the different zones in the modeland their corresponding acoustic impedance value - this is the case where a well has been drilled through the salt bodyas shown in Figure 5c. Similar to the previous example, the post-stack seismic data (Figure 5b) is modelled using the G operator with a 8Hz Ricker wavelet.Figure 7 shows the estimated model for three different solvers: LSQR (used for the least-squares regularized functional),L-ADMM and PD (used the TV regularized function). Once again, the inversion of noisy data benefits from use ofthe isotropic TV norm which regularizes the solution and produces sharp discontinuities at the edges of the model.Nevertheless, a common feature of the three inverted models is represented by the inaccurate estimate of the acousticimpedance of the salt body which is underestimated by all of the tested algorithms. This is not surprising as thebackground model is fairly far from the true solution, especially towards the edges of the salt body, and the frequencycontent of the seismic data does not allow to recover such a gap even when a strong regularization such as the isotropicTV norm is included in the functional to optimize. 8oint Inversion-Segmentation for Seismic Interpretation A PREPRINT x[m] t [ s ] (a) True x[m] (b)
Segmentation from PD x[m] (c)
Detected Horizons (d)
Probs.
Figure 5: Same as Figure 3 with the Primal-Dual model used as input for segmentation. x[m] t [ s ] (a) AI[kg/(m s)] x[m] (b) Seismic (c)
AI trace
Figure 6: a) Acoustic impedance model, b) noisy seismic data, and c) AI profile at well location (shown by the dashedblack line in panel a)Despite of the observed inaccuracies in the retrieved acoustic impedance model, segmentation and horizon trackingstill produce satisfactory results when applied to the Primal-Dual model (Figure 7). Note that most of the zones in theoriginal model (Figure 8a) have been segmented precisely (Figure 8b) with only the shallowest and two deepest zonesbeing quite noisy due to the weak impedance contrast affected by noise in the data (Figure 6b). A similar conclusion canbe drawn for the outcome of the horizon tracking step (Figure 8c) where most of the horizons overlay almost perfectlywith their true counterpart apart from those related to the poorly recovered zones.The joint inversion-segmentation algorithm proposed in the previous section is now employed to evaluate whetherfurther constraining repeated steps of inversion with previously segmented zones could further improve the results inFigures 7 and 8. The retrieved model in Figure 9a after 4 outer iterations of the joint algorithm does indeed more closelyresemble the true model. This is especially the case for the salt body, whose absolute value is not underestimated in thiscase (Figure 10a), as well as the overall continuity of the different layers. This is the consequence of the fact that thefirst step of segmentation is used as soft constraint to the second step of inversion (last term in equation 2a), whoseoutput is in turn used to drive the second step of segmentation (first term in equation 2c) and so on and so forth. Theoutcome of the last segmentation step (Figure 9b) is also clearly more accurate than the one produced by directly usingthe Primal-Dual model: we observe a better continuity of the top interface of the shallowest zone as well as closer9oint Inversion-Segmentation for Seismic Interpretation A PREPRINT x[m] t [ s ] (a) L2 (PSNR=50.20) x[m] (b) L-ADMM (PSNR=50.96) x[m] (c)
PD (PSNR=51.81)
Figure 7: Inverted models via a) Least-squares regularized inversion, b) Anisotropic TV-regularized inversion Split-Bregman solver, and c) Isotropic TV-regularized inversion Primal-Dual solver. d) True, background, and inverted AIprofiles at well location. x[m] t [ s ] (a) True x[m] (b)
Segmentation from PD x[m] (c)
Detected Horizons (d)
Probs.
Figure 8: a) True and b) estimated segmentation from Primal-Dual model, respectively. c) Tracked horizons overlaid tothe noise free seismic dataset (thin lines: true horizons, thick lines: tracked horizons. d) Probabilities for each class atwell location (i.e., a single column of V at the index in the middle of the model.resemblance to the true classes in the deeper parts of the model. Similarly, the tracked horizons are of overall higherquality especially on the left of the salt body and in the deeper part of the model.A second scenario is now evaluated. The acoustic impedance model in Figure 6a is modified by adding small scalefluctuations to the macro model composed of a limited number of acoustic impedance values. Such small scale variationsare clearly visible in the log data (Figure 10b) and further complicate the overall seismic response (Figure 11e). Thedefinition of our classes c is therefore made by blocking the acoustic impedance log. We clearly do not expect ourinversion to be able to recover such features given that they are outside of the bandwidth of the signal: on the other hand,this examples serves the purpose of testing the sensitivity of our algorithm to more complex and realistic models of thesubsurface. Panels c, e, and f in Figure 11 show that the different steps of our joint inversion are robust to small scalevariations in the model and the results are overall very similar to those of the ideal case with a well defined macro model.Whilst all the horizons are successfully tracked in the segmented model, some spurious events are also reconstructed aspart of the deeper horizon. This is the direct consequence of small blobs in the deeper layer of the segmented model(11e) which lead our tracking algorithm to identify those lines as part of the same horizon group as the deepest horizonbecause they border same pair of classes as the larger horizon does.Finally, we investigate the robustness of our algorithm with respect to partial knowledge of the different classes. Asshown in Figure 12 the vertical well is now assumed to penetrate only some of the layers in the model: as a consequence,the classes vector c only contains seven elements to which we add an eighth based on the assumption that we are awareof the presence of the salt body and have a good estimate of its acoustic impedance. Whilst our algorithm fails to trackthe two deepest interfaces, both the inversion and horizon tracking of the shallow subsurface is still satisfactory and notaffected by the lack of knowledge of the deeper part of the model. On the other hand, as shown in Figures 12a and10oint Inversion-Segmentation for Seismic Interpretation A PREPRINT x[m] t [ s ] (a) Joint (PSNR=53.43) x[m] (b)
Segmentation from Joint x[m] (c)
Detected Horizons (d)
Probs.
Figure 9: a) Model and b) segmentation from joint inversion scheme. c) Tracked horizons overlaid to the noise freeseismic dataset (thin lines: true horizons, thick lines: tracked horizons. d) Probabilities for each class at well location(i.e., a single column of V at the index in the middle of the model.Figure 10: Acoustic impedance estimates for the different inversion algorithms at well location for Hess dataset. a)Blocky model in Figure 6a, b) Model with small scale perturbations in Figure 11a, and c) Model with short well inFigure 12a.Figures 10c, the acoustic impedance estimates in the deeper section are underestimated as a consequence of the fact thatour joint algorithm tends to drive acoustic impedance values at each location closer to those of the class that has beenselected in the previous iteration. Volve is an oil field located in the central part of the North Sea, five kilometres north of the Sleipner Øst field. Volveproduced oil from sandstone of Middle Jurassic age in the Hugin Formation, with the main reservoir located at a depthof approximately 2,700-3,100 metres. The field was shut down in 2016, with the facility removed in 2018, and allhistorical subsurface and production data made available by Equinor in June 2018.In this section, our joint inversion and segmentation algorithm is applied to a 2d section of the PSDM full stack datasetfrom the ST10010ZC11 survey (Figure 13a). Since interpretation has been carried out in time, the time version of thisdataset is used in our example. The section is extracted along the NO/15-9 19 BT2 well and further extended to theEast of the well as shown on top of a time map of the BCU time surface (Figure 13b). Figure 13c shows the acousticimpedance log converted into two-way traveltime (TWT) using the available checkshot profile. The 3 key interpretedhorizons are overlain on the seismic section in Figure 13a and their corresponding well markers are shown in Figure13c. 11oint Inversion-Segmentation for Seismic Interpretation A PREPRINT t [ s ] (a) True (b)
PD (PSNR=51.81) (c)
Joint (PSNR=53.60) x[m] t [ s ] (d) Seismic x[m] (e)
Segmentation from Joint x[m] (f)
Detected Horizons (g)
Probs.
Figure 11: Results for model with small scale variations. a) True model, b) model from Primal-dual inversion, c) modelfrom joint inversion scheme, d) noisy data, e) segmentation from the joint inversion scheme, f) tracked horizons overlainon the noise-free seismic dataset (note that noisy features are due to small scale AI perturbations in the model), and g)probabilities for each class at well location. x[m] t [ s ] (a) Joint (PSNR=49.55) x[m] (b)
Segmentation from Joint x[m] (c)
Detected Horizons (d)
Probs.
Figure 12: Results for model with short well. a) Inverted model from joint inversion scheme, b) segmentation from thejoint inversion scheme, and f) tracked horizons overlain on the noise-free seismic dataset.12oint Inversion-Segmentation for Seismic Interpretation A PREPRINT
To begin with, the acoustic impedance profile is used to identify classes for segmentation step. Three different zones arechosen, separated by the Shetland and Viking tops and the histograms of their acoustic impedance values plotted inFigure 14. However, given the large overlap in the histograms from the overburden class and the Viking class, they areaggregated into a single class. A second class is defined to include most of the values in the Shetland formation and twoother classes are defined to capture the low and high acoustic impedance values observed in the well log. Ultimately,4 different classes ( c = [4000 , , , ) are used as input to our segmentation algorithm. Note that ourability to successfully separate different formations is highly dependant on the choice of these classes. To accuratelyfine tune such a step, it may be useful to perform a first independent inversion step and apply a pixel-wise segmentationbased on the closeness of each pixel value to the classes values. This is equivalent to solving the segmentation problemin equation 2c without the TV-regularization term ( β = 0 ) leading to a very noisy segmentation result. Nevertheless, byclosely inspecting this segmentation the classes boundaries can be further optimized prior to running the entire jointinversion and segmentation scheme.To be able to invert the seismic data for an acoustic impedance model, a background model is built from the root-mean-square (RMS) velocity model. RMS velocities are first converted into interval velocities and subsequently calibratedwith the acoustic impedance log of the NO/15-9 19 BT2 well. More specifically, the interval velocity model is extractedalong the well trajectory and used to ‘predict’ a smoothed version of the acoustic log: prediction is achieved via linearregression and the regression coefficients are further used to convert the entire velocity model into a background acousticimpedance model (Figure 15a). Note that other approaches to construct the low-frequency model can be equivalentlyused in cases where a larger coverage of vertical wells with the required set of logs is available. Estimated acousticimpedance models for three different optimization problems are shown in the other panels of Figure 15: from left toright, spatially regularized least-squares inversion, TV regularized inversion, and our joint inversion-segmentation. Asalready observed in the other examples, TV regularized inversion improves on the blockiness of the model compared toleast-squares inversion and produces sharper transitions between different formations; this is especially the case forareas with low and high acoustic impedance values which are generally under predicted due to the lack of low frequencyinformation in the data. After two outer iterations of the algorithm in equation 2, this behaviour is even more visible asa consequence of the second regularization term in equation 1, which drives acoustic impedance values closer to that ofthe class they in which they belong to. Moreover, whilst the other two inversion results show some vertical striping inan area of poorer data quality (white arrows), the joint inversion manages to improve the lateral continuity of the modelwithout compromising on the sharpness of the vertical transitions. Finally, an overall good match is observed with theacoustic impedance log along the well trajectory as shown in Figure 15e.Figure 17 displays the segmented model at the end of the joint inversion scheme alongside with the probabilities of eachclass along the well trajectory. White lines represent the different horizons that have been extracted from the differentclass probabilities. Due to the complexity of the model and the fact that some pairs of classes repeat at different depths,the horizon tracking algorithm is run in this case only considering the class above during the labelling and combinationsteps. The tracked horizons are also shown in Figure 17c alongside the manually interpreted horizons (thin blacklines) on top of the input seismic data. A number of interesting observations can be made with respect to the ability ofour algorithm to interpret the key horizons in this data: first, the Top Shetland formation is successfully tracked andthe resulting horizon is very similar to the manually interpreted one. On the other hand, our interpretation of Top Tydiverges from the manual interpretation towards the right side of the model; a similar behaviour is also observed for thehorizon above (red line in Figure 17c). By looking at the inverted acoustic impedance (Figure 18b), we can observe athinning of the Ty formation which is consistent with the tracked horizons. Whilst the quality of the seismic data in thisarea does not allow to determine whether our interpreted horizons are correct, this result highlights the direct connectionbetween the tracking algorithm and the inversion step. Future work will investigate this behaviour by both invertingthe entire 3d dataset and looking at its consistency in the perpendicular direction as well as looking at the pre-stackdata in this area. Finally, our algorithm is only able to track part of the BCU horizon, failing to do so when the thinlow acoustic impedance layer is eroded out on the right side of the model. Overall, the results in this section confirmthe validity of our algorithm and its usefulness in jointly solving a number of tasks and keeping consistency amongthe different results: in other words, whilst the acoustic impedance model alone may be useful to condition facies andproperty modelling, the segmentation model may complement it in defining different areas of influence which may beconditioned differently. Discussion
The computational cost associated with the Primal-Dual algorithm used to optimize equations 2a and 2c is highlydependant on that of the proximal operators ( prox f and prox g ) involved in each iteration. As already shown in equation13oint Inversion-Segmentation for Seismic Interpretation A PREPRINT
Figure 13: a) Seismic section along the NO/15-9 19 BT2 well with Top Ty (brown line), Top Shetland (blue line) andBCU (red line) horizons. b) BCU time horizon with well trajectory and fence along which seismic is extracted anddisplayed in the previous panel. c) Acoustic impedance well log and well markers for the 3 formations of interest. Alldata displayed here is taken from the official Volve dataset.
AI distribution
OverShetlandViking
Figure 14: Histogram of the acoustic impedance values at NO/15-9 19 BT2 well divided in zones as defined in Figure13c. Solid black vertical lines refer to the classes chosen for the segmentation algorithm whilst dashed black verticallines represent the separation between the classes.9, the proximal operator of the L , requires only point-wise evaluations for each element of the input vector; this istherefore very fast to compute even for very large two or three dimensional models.On the other hand, the evaluation of the proximal operator for the quadratic functional || d − Gm || requires thesolution of a regularized least-squares inverse problem. Apart from some specific cases involving orthogonal operators,where an analytical solution can be identified (e.g., CT scan reconstruction - [12]), the evaluation of such an operator isgenerally very expensive and requires the use of iterative solvers (e.g., CGLS, LSQR). However, as each iteration of thePrimal-Dual algorithm represents only a step towards the saddle point of the corresponding Primal-Dual functional, anapproximated solution of the problem obtained with an early stopping after a limited number of iterations generallysuffices. Moreover, a warm start is generally employed, meaning that the solution of the proximal operator fromits previous evaluation is used as starting guess of its current evaluation. This often provides a very large speedimprovement over solving the problem from scratch each time [14].Finally the proximal operator of the unit Simplex used within the segmentation process (equation 2c) also requires onlypoint-wise evaluations, where each column of the matrix V is treated independently and projected over the intersectionof a hyperplane and a box. This is implemented as the evaluation of a projection over a box: P C = P Box [0 , + inf] ( x − µ ∗ ) (11)where µ is the solution of the following scalar equation f ( µ ) = T P Box [0 , + inf] ( x − µ ) − (12)14oint Inversion-Segmentation for Seismic Interpretation A PREPRINT T W T [ s ] (a) Background (b) Least-squares T W T [ s ] (c) PD (d) Joint T W T [ s ] (e)Along well AIlog(t)BackL2PDJoint
Figure 15: Acoustic impedance models for the Volve dataset: a) background, b) regularized least-squares, c) primal-dual,and d) joint inversion. e) AI values at well location for the different inversions.
350 400 450 500 550 T W T [ s ] (a) Least-Squared
350 400 450 500 550 (b)
Primal-Dual
350 400 450 500 550 (c)
Joint
Figure 16: Zoomed sections of acoustic impedance inversion for a) least-squares solver, b) Primal-dual solver and c)Joint inversion. 15oint Inversion-Segmentation for Seismic Interpretation A PREPRINT T W T [ s ] (a) Segmentation (b)Probs. T W T [ s ] (c) Detected Horizons Figure 17: a) Segmentation model with tracked horizons (white lines), b) seismic data with tracked horizons, and c)probabilities for each class at well locationwhich can be solved via bisection. The cost of this operation is dependant on the speed of convergence of the bisectionalgorithm and may require tens of evaluations of the function f . Note that, although this has a higher cost comparedto the proximal operator of the L , norm, its implementation can be made very fast when implemented on modern,massively parallel hardware (i.e., GPUs). Extension to pre-stack data
The algorithm presented in this work can be easily adapted to the case of pre-stack seismic data with the possiblyadded benefit of better separation between the different classes chosen as input to the segmentation step. As far as theinversion step is concerned, the model vector m can simply be defined as a vector of size N m N x N z × where N m is the number of model parameters we wish to invert for (generally 2 or 3). Similarly, we can define the data vector d as a vector of size N θ N x N z × where N θ is the number of angles. Finally, the modelling operator is redefined as G = WMD to include a mixing matrix M which contains the weighting coefficients of the different elastic reflectivitiesused to model the amplitude variation with offset (AVO) seismic response - e.g., Aki-Richards. Similar to the matrix V , we can also consider a matrix M of size N x N z × N m (where m = V ec ( M ) is its vectorized version obtained byconcatenating the matrix’s columns) and replace the TV regularization term in equation 1 by the sum of N m TV normsacting on each column of M . The proximal operator of this term can therefore be computed in the same way as that of16oint Inversion-Segmentation for Seismic Interpretation A PREPRINT T W T [ s ] (a) (b)
200 300 400 500 600 700 T W T [ s ] (c)
200 300 400 500 600 700 (d)
Figure 18: Zoomed sections of (a-c) seismic data and (b-d) acoustic impedance inversion with interpreted horizons (top:original interpretation, bottom: our interpretation). V (equation 9). Alternatively, as described in [18, 19], the different model parameters can be first decoupled from eachother, followed by N m independent steps of inversion that are identical to the one employed in the post-stack scenario.Finally, for the segmentation step we can define a matrix C of size N c × N m whose rows contain the set of elasticparameters associated to each class and replace the ( m i − c j ) in equation 1 by the Euclidean norm of the differencebetween each row of M and || M i − C j || . Conclusions
In this work we have presented an end-to-end approach to assisted seismic interpretation which provides alongside a setof interpreted horizons also a property model (or a set of properties models for the elastic case) and a segmented modelof the subsurface. Acoustic (or elastic) properties are not only inverted for the entire model from seismic data, but theyare also used as input to define an easy to understand link between well log responses and the expected horizons to beextracted in the final step of our workflow. This in turn reconciles classical principles of seismic interpretation to theproposed optimization problem, making the method easier to explain to and be used by skilled interpreters. Our workremarks on the importance of enforcing blockiness in the inversion of seismic data, which is obtained via a combinationof regularization (TV regularization), solver selection (Primal-Dual) and additional constraints (segmentation model):our estimates present sharper jumps at model discontinuities than their least-squares counterpart. This not only helpsthe subsequent segmentation step, but it has strong implications when these models are used as hard or soft constraintsin downstream processes such as facies or property modelling. Finally, whilst an extension of our inversion algorithm topre-stack seismic data has been discussed, we foresee its applicability to other linear (and nonlinear) inverse problemsin geophysics, such as least-squares migration and waveform inversion.
References [1] D. A. Herron.
Seismic Intepretation . Society of Exploration Geophysics, 2011.17oint Inversion-Segmentation for Seismic Interpretation A PREPRINT [2] A. U. Waldeland, A. C. Jensen, L.-J. Gelius, and A.H.S. Solberg. Convolutional neural networks for automatedseismic interpretation.
The Leading Edge , 37(4):529–537, 2018.[3] D. Haibin, L. Truelove, C. Li, and A. Abubakar. Accelerating seismic fault and stratigraphy interpretation withdeep cnns: A case study of the taranaki basin.
The Leading Edge , 39(10), 2029.[4] Y. Shi, X. Wu, and S. Fomel. Saltseg: Automatic 3d salt segmentation using a deep convolutional neural network.
Interpretation , 7(3):SE113–SE1221, 2019.[5] Y. Shi, X. Wu, and S. Fomel. Waveform embedding: Automatic horizon picking with unsupervised deep learning.
Interpretation , 85(4):1–48, 2020.[6] D. Hampson, B. Russell, and B. Bankhead. Simultaneous inversion of pre-stack seismic data.
SEG TechnicalProgram Expanded Abstracts , pages 1633–1637, 2005.[7] S. Osher, M. Burger, G. Goldfarb, J. Xu, and W. Yin. An iterative regularization method for total variation-basedimage restoration.
Multiscale Model. Simulation , 4:460–489, 2005.[8] T. Goldstein and S. Osher. The split bregman method for l1-regularized problems.
SIAM Journal on ScientificComputing , 2(2):323–343, 2009.[9] A. Gholami. Nonlinear multichannel impedance inversion by total-variation regularization.
Geophysics , 80:R217–R224, 2015.[10] D. Wang, J. Gao, and H. Zhou. Data-driven multichannel seismic impedance inversion with anisotropic totalvariation regularization.
Journal of Inverse and Ill-Posed Problems , 2018(2):229–241, 2019.[11] O. Kolbjornsen, O. Evensen, A. K. Nilsen, and J. E. Lie. Digital superresolution in seismic avo inversion.
TheLeading Edge , 39:791–799, 2019.[12] V. Corona, M. Benning, M. Ehrhardt and. L. Gladden, R. Mair, A. Reci, A. Sederman, S. Reichelt, and C.-B. Schonlieb. Enhancing joint reconstruction and segmentation with non-convex bregman iterationt.
InverseProblems , 2019.[13] M. Benning and M. Burger. Ground states and singular vectors of convex variational regularization methods.
Methods Appl Anal. , 20:295–334, 2013.[14] N. Parikh.
Proximal Algorithms . Foundations and Trends in Optimization, 2013.[15] A. Chambolle and T. Pock. A first-order primal-dual algorithm for convex problems with applications to imaging.
Journal of Mathematical Imaging and Vision , 40:120–145, 2011.[16] M. Ravasi and I. Vasconcelos. Pylops - a linear-operator python library for scalable algebra and optimization.
SoftwareX , 11:100361, 2020.[17] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed optimization and statistical learning via thealternating direction method of multipliers.
Foundations and Trends in Machine Learning , 3(1):1–122, 2011.[18] E. Causse, M. Riede, A.J. Van Wijngaarden, A. Buland, J.F. Dutzer, and R. Fillon. Amplitude analysis with anoptimal model-based linear avo approximation: Part i-theory.
Geophysics , 72(3):C59–C69, 2007.[19] M. Ravasi, M. Alerini, S. Maultzsch, and A. Ghaderi. Band-limited optavo, seismic inversion the other way round.79th EAGE Conference and Exhibition