History-based Interpretation of Finite Element Simulations of Seismic Wave Fields

1 Introduction

Most work on qualitative reasoning about physical sys­ tems is devoted to technical systems consisting of a fixed set of components that interact via given connections. Examples given in [Weld and de Kleer, 1990] are elec­ tronic circuits, water tanks , and gear systems. In con­ trast to this, we will consider natural systems where ob­ jects are dynamic in position, direction, shape, and num­ ber. The FROB system [Forbus, 1984] simulates spring­ ing balls changing their positions and directions, but keeping their shapes. Furthermore, we don 't obtain new balls. Flowing liquids [Hayes, 1985a] are different : They easily divide, merge , and change their shapes. In order to capture those interactions between liquids, Hayes de­ veloped the concept of a history, i.e. a coherent piece of space-time. Histories provide an adequate means to de­ scribe the behaviour of dynamic objects such as flames, waves, clusters, clouds , which can all be deformed , di­ vided and merged.

In this paper, we will consider a concrete task requir­ ing history-based reasoning about physical phenomena. We consider the propagation of seismic shock waves in the underground [Lavergne, 1986] . Seismic waves are used by geophysicists to explore the structure of the underground. They are usually launched by an initial vibration on the surface. The resulting spheric shock wave is then propagating downwards as shown in the first snapshot of figure 1. When it hits an interface be­ tween two geological layers this causes a reflected and a

·This paper is based on work performed during the post­ doctoral stay of the first author at the lnstitut Fran~ais du Petrole.

Bertrand Braunschweig

lnstitut Fran<_;ais du Petrole 1 et 4, av. de Bois-Preau, BP 311 92506 Rueil-Malmaison France

Bertrand.Braunschweig@ifp.fr

transmitted wave. The reflected wave returns to the sur­ face and leaves an observable front in the seismograms measured by the geophysicists.

In order to interpret seismograms, the geophysicists incrementally construct a model of the underground based on hypotheses of the histories of the returning waves. Above, we considered a wave that was reflected by the first interface. Further interfaces lead to further direct reflections . Additionally, a seismogram can show multiply reflected fronts, diffractions which are obtained due to corners and many other disturbing fronts . Geo­ physicists pick out direct reflections using some heuristic approach and use them to construct a model of the un­ derground (based on numerical optimization procedures or further ad-hoc rules) .

Newer work on numerical simulators based on finite elements allows a very precise simulation of the wave propagation in complex models of the underground . The snapshot sequence in figure 1 has been produced by such a simulator lAnne and Brae, 1994]. The simulations en­ able a verification of the geological model. Divergences between observed and simulated seismograms might help to correct the model. To detect them, we have to com­ pare fronts having the same history (e.g. two direct re­ flect ions; two diffractions etc.). Unfortunately, numeri­ cal simulators based on finite elements do not keep track of the history of waves . They produce a series of images showing the waves, but they lack a representation of the wave objects, their causes, and their histories. When ex­ amining a front of a seismogram, we want to know the obstacles and the types of phenomena that produced it.

In this paper, we show how to interpret the images pro­ duced by the numerical simulator and how to establish a causal relation between seismic events , waves, and obsta­ cles in the underground . Our goal is to detect Hayes-like histories of waves in snapshot images. Although the pa­ per is restricted to 2D-models of the underground , its concepts can be generalized to the 3D-case.

The paper is divided into two main sections . Section 2 presents the representation of fields (sec. 2.1), as well as the vocabulary for describing wave histories (sec. 2.2) . The interpretation is done in several steps developed in section 3. We first decompose the underground into lay­ ers and interfaces (sec. 3.1) . Then we show how to detect wave fronts in a single snapshot (sec. 3.2). Tracking a front from one snapshot to the other is considered in sec­ tion 3.3. In section 3.4, we discuss how to detect new objects and their causes.

2 Multiple representations

2.1 Fields

In order to describe complex phenomena (e.g. liquids, waves, flames etc.), physicists use parameter fields . A field is the distribution of a physical parameter in the given space. For shock waves, we consider a velocity field, giving the velocity of a wave at a certain point, and the field of the amplitudes of the waves (i .e. the snapshots in figure 1.). A physical law captures a rela­ tionship between parameter fields, which is valid at each point. In general, such a law is a differential equation (e.g. the wave equation). Its solution describes the tem­ poral development of a field. Some of the fields such as the velocity field are static (stationary), whereas the amplitude field is changing in time (non-stationary). We restrict our discussion to a single static and a single dy­ namic field.

A well-suited technique for simulating changes of com­ plex and arbitrary fields is the finite element method. A numerical simulator based on this technique is supplied with the initial parameter fields and then produces a se­ ries of snapshots showing the fields at selected instances of time. The finite element method can handle complex shapes because it uses a kind of an analogical representation of fields: '

• It is based on a grid (P, N) where P is a set of se­ lected points and N ~ P x P is a symmetric neigh­ bom·hood relation.
• It describes the spatial distribution of a parameter using a mapping f : p ....... n of the points to the real numbers.
• It is specified extensionally (e.g. a matrix of floating point numbers.)

Seismic simulators normally use regular grids obtained by rows and columns. They are characterized by a start­ ing point s := (s 1 , s2), a unit distance Ll, the number n of columns and the number m of rows. The set of points is then given by

Two points are neighbours if they have successive posi­ tions in the same row or column.

2.2 Histories of dynamic objects

Fields don't represent objects explicitly. They just show certain patterns of activity that are reproduced in the next instants. For example, figure 1 shows wave fronts that are propagating, hitting interfaces, and generating new waves. In order to describe these phenomena, we need an ontology for dynamic objects in fields.

Our discussion is based on a given (continuous) space s, for example the two-dimensional space defined by n 2 , and a linear (continuous) time defined by T := R. Dy­ namic objects such as waves evolve in time and occupy a region at each timet. This region is a subset of S x {t}. If we consider different time points the occupied region of an object can change. We require that these changes are local. If we put the regions of an object at different times together, we obtain a subset of S x T . This subset must be a 'connected piece of space-time', i.e. a history as defined in [Hayes , 1985b].

The region occupied by an object can change in a con­ tinuous or discontinuous way. For example, the initial wave front in figure 1 is split into two parts when hit­ ting the interface. Changes are caused by the dynamic and the static field. In the case of seismic waves, discon­ tinuities caused by the static field can be characterized precisely: If the velocities of a zone are changing contin­ uously, a (convex) wave propagating through this zone will change continuously. Discontinuities in the velocity field however cause discontinuities in wave propagation.

In order to keep track of causes for discontinuous changes, we divide the static field into regions and ob­ stacles. A region is a (maximal) coherent subset of S that does not contain discontinuities in the static field. The regions of the seismic velocity field are also called geological layers. An obstacle separates one, two, or sev­ eral regions. It is a (maximal) coherent subset of S of discontinuity points in the static field. Its dimension is inferior to that of S. The geological model of figure 1 is composed of three 2D-regions of constant velocities, which are separated by three 1D-obstacles called inter­ faces. The interfaces are linked by a corner which is a OD-obstacle. Regions, interfaces , and corners consti­ tute a place vocabulary in the sense of [Hayes , 1985b; Forbus, 1994].

Thus, we divided the static field into regions where motion of dynamic objects is continuous and into obsta­ cles which disturb motion in a discontinuous way. We now use this topological structure to divide the global histories into episodes of dynamic objects and incident events linking these episodes. We require that a dynamic object is contained into a single region. If it reaches an interface then the continuation of its history on the other side of the interface is considered to be a new object, namely the transmission of the incident object. A dy­ namic object is an episode of the global history that is contained in the static history of a certain region.

An incident is the event when a dynamic objects hits an obstacle. It is in fact the intersection of the history of the dynamic object with the history of the obstacle. The incident is caused by the incident object and causes itself new objects in the regions surrounding the obstacle. The global history is branching at the incidents as illustrated in figure 2. An incident is the start of the histories of the waves it causes and it links them with the history of the incident wave

Thus, we have structured histories in the dynamic fields into dynamic objects, incidents, and their causal relation. In fact, we have adapted the basic concepts of naive physics [Hayes, 1985b] to physical fields and now have a vocabulary for interpreting the simulated fields.

3 Interpreting snapshots

3.1 Detecting static histories

First we show how to decompose the static velocity field into geological layers and interfaces. Thus, we obtain the regions where to look for waves and the obstacles where to look for incidents.

Let (P, N) be a grid and f : P -+ R be a field. We define regions as follows: let C ~ N be a symmetric cri­ teria that specifies whether two neighbour points belong to the same region. We consider the reflexive transitive closure of c and denote it by c·. c· is the smallest superset of C that is reflexive and transitive. Since C is symmetric, c· is an equivalence relation . The regions are obtained as the equivalence classes of c·. The C­ region of a point p E P is defined as the equivalence class containing p:

To define regions in the velocity field v of seismic waves, we link two neighbour points if there is no discontinuity between them. Since grids have a fixed resolution, we use a threshold c to operationalize this criteria. The velocity difference of two points must be smaller than c:

The geological layer of point p is then the equivalence class Rv(p).

Next we define interfaces separating two C-regions R1 and R2 . An interface isjust a set ofneighbourhood links (p, q) E N that do not satisfy the given criteria C and that connect a point in R1 with a point in R2 :

The set I(R1 , R2) is called C-interface between R1, R2 iff I(R1 , R2) is not empty. The geological layers are sepa­ rated by V -interfaces. For the sake of brevity, we neither discuss corners, nor the case that the interface between two regions is interrupted by a third region.

3.2 Detecting objects in a snapshot

In the next sections, we consider a sequence of snapshots a1 , a2 , a3 , ... showing the amplitude field at increasing time points t 1 , t2, t 3 , . . .. We proceed in three steps in order to detect histories of wave objects. First, we iden­ tify wave objects in a single snapshot. Then, we link the possible interpretations of succeeding snapshots. Af­ ter that, we show how to detect histories of new objects caused by incidents at interfaces.

Wave fronts as shown in figure 1 consist of a small number of oscillations. In a snapshot, they appear as thin regions of negative or positive amplitudes, which can clearly be distinguished from the background having zero amplitude. To capture this phenomena formally, we divide the set P of points into three classes: positive, negative, and zero ones. Since there are small distortions in the simulated field, we use a 5 > 0 to define the zero class. Let a; be the amplitude field of the i-th snapshot:

Coherent regions of either positive, negative, or zero points are obtained by linking two neighbour points of the same class. Since we want to reconstruct wave fronts and each wave front is limited to a single geological layer we cut the regions in the amplitude field at the geological interfaces by using the criteria V.

This criteria A; gives rise to A;-regions and A;-interfaces decomposing the geological layers at time i; into positive, negative, and zero regions.

Characterizing wave fronts by their positive and neg­ ative regions is problematic because we can obtain in­ tersections between several fronts ( cf. figure 4). Due to this, the regions of several waves can merge. A more elegant and robust characterization is obtained by the zero regions that are enclosed by a wave. If a wave front f in the layer L is enclosing a zero region A then we obtain the front as the interface between the regions A and L-A. In case of interruptions, a front encloses sev­ eral zero regions. Front fi in figure 4 encloses the zero regions A, D, whereas front h encloses the regions D, C. fi also encloses the gap between the fronts A and D which is caused by f2.

In order to close this gap, we will enlarge the zero regions of a layer L until the complete layer is covered. A single growing step adds the neighbouring layer points to a region X~ P:

We iterate this step until the layer L is completely cov­ ered. Let Z;(L) be the set of the zero regions contained in the layer L at time i;. The i-th cover factor of L is the smallest number k s.t.

If Z is a zero region in layer L at time i; and k is the cover factor of L at this time, then the enlargement of Z is defined as

The enlargement enable us to define a wave front using the zero regions enclosed by it . Suppose that w is a wave in layer L and that the front of w at time t; encloses 1 the zero regions Z1 , ... , Zm of L. Then

1 Our characterization is based on the assumption that

1. Z;(w) := {Z1 , .. . ,Zm} is called i-the characteriza­ tion of w.
2. from;(w) := Z1 u ... UZm is called the i-th enclosed region of w.
3. to;(w) := L- from,(w) is called the i-th propagation region of w.

The interface between from; ( w) and to; ( w) represents the wave front. For example, the front fi is obtained as the interface between AU band B U CUE, whereas h is the interface between CUb and B U AU E.

If from;(w) is empty then the wave w has not yet ap­ peared in snapshot i. If to;(w) is empty w has gone.

3.3 Tracking existing histories

Dynamic objects are steadily changing form and posi­ tion. The time steps between two succeeding snapshots are usually too large to track these changes locally ( cf. figure 1). In this section, we show how to track the his­ tory of an object even in presence of gaps.

For a given set W of waves and the set Z of zero re­ gions of snapshot i + 1, the task is to identify the i + 1-th characterizations of the waves in W by using their i-th characterizations. We require that the new characteri­ zations satisfy some additional criteria:

1. The new fronts should be in the propagation regions of the old fronts.
2. Each observed frontier between two zero regions should be explained by a wave.
3. The new fronts should be as close as possible to the old fronts.

The first constraint is based on the hypothesis that all points are propagating to the same side of the front. It implies that the enclosed region of a wave w is grow­ ing from one instant to the other, which is expressed by from; ( w) ~ fromi+l ( w). If a zero region Z E Z of snapshot i + 1 is overlapping with from; ( w) it is also overlapping with from;+ 1 (w) and therefore an element of the i + 1-th characterization of w:

The second constraint requires the detection of fron­ tiers between two zero regions in a snapshot. In figure 4, D has a frontier with A and C, but not with B and E. If we let grow all regions then the borders of D will overlap with that of A, B, C, E. However, the overlap with B and E is also covered by the borders of A and C. This leads to the following definition: There is a frontier between two zero regions X and Y in Z iff

A frontier between X and Y is explained by a wave w iff w is enclosing exactly one of these zero regions, i.e. I {X, Y} n Z;+l(w) I= 1.

wave fronts always enclose a region. If a front ends inside a layer we cau further divide the zero regions by lengthening this end.

These two constraints already reduce the number of possible interpretations. Consider figure 4 showing two waves WI and w2 at two snapshots 1 and 2. Their wave fronts in snapshot 1 are indicated by dotted lines. The 1-st enclosed region of wave WI overlaps with the zero re­ gions A and D, whereas the 1-st enclosed region of wave w2 overlaps with C and D. Due to our first constraint, we get:

In order to explain the frontier between C and D, the region C cannot be an element of Z2(wi) · Similarly, A cannot be an element of Z2(w2 ). What about Band E? Their frontiers to the regions A and C are explained if they are not enclosed by wave WI and w2. However, the frontiers of B are also explained if B is enclosed by WI and w2 . A similar argument holds for E. Thus, there are four characterizations explaining all frontiers:

In order to reduce these ambiguities, we require that wave fronts are as close as possible to the old fronts. This can be achieved by minimizing the characterizations of waves. Let Zl+I : W --. 2 2 and zp+I : W --. 2 2 be two characterizations. Zl+I is preferred to Zl+I iff

for each w E W. We don't accept a characterization Zf+I if there exists a characterization Z;*+l that is pre­ ferred and different to z:+I and that satisfies the two constraints introduced above.

This eliminates the characterizations 2, 3, 4 of our ex­ ample. Hence, the wave WI is characterized by {A, D} and the wave w2 is characterized by { C, D} . Our three principles are sufficient to track histories if the gaps be­ tween two succeeding fronts of a wave are not too large.

3.4 Detecting new histories

Intersections between different histories can be the start of new histories. In the case of seismic waves, the in­ tersection of two wave histories is without interaction, whereas the intersection between a wave history and the static history of an obstacle bears new wave histories in the regions surrounding the obstacle. In this section, we show how to detect incidents and their resulting waves. For the sake. of shortness, we consider only incidents at interfaces causing reflections and transmissions.

First, we discuss how to detect and characterize inci­ dents to interfaces. Let w be a wave in layer LI and I be the (coherent) interface Iv(L1 , L2 ) between £ 1 and a neighbouring layer £2. The event of an incident of w to I is denoted by e(w, I).

To detect an incident e( w, I) in snapshot i, we are checking whether w hits I in this snapshot. A wave hits an interface iff the region behind it, as well as the region in front of it are touching this interface. A region X s; P touches an interface I s; N if some points of X are linked to other points by I . In this case, the following set of edges is non-empty:

We can now characterize incidents by the points in to; ( w) and from;(w) that are touching I:

The set from; ( e( w, I)) is called the i- th enclosed interface of the incident e(w, I), whereas to;(e(w, I)) is called the i-th propagation interface of the incident e( w, I) .

If the i-th enclosed interface is empty then the incident e( w, I) has not yet occurred in snapshot i . If the i­ th propagation interface is empty then the incident is completed in snapshot i. If neither the i-th propagation interface, nor the i-th enclosed interface are empty then the incident occurs in snapshot i.

If an incident e := e( w, I) occurs in a snapshot i then it causes a reflected wave 1·( e) in the layer LI of the incident wave and a transmitted wave t( e) in the layer £ 2 on the other side of the interface. These waves exist if and only if the incident occurs:

The enclosed region of the reflected wave is touching the interface I exactly at the enclosed interface of the inci­ dent . Furthermore, the enclosed region of the reflected wave is included in the enclosed region of the incident wave:

The transmitted wave is touching the interface I at least at the enclosed interface of the incident:

Vle don't get the inverse inclusion because the transmit­ ted wave can be fa.'3ter than the incident wave.

VVe are thus able to detect reflections and transmis­ sions of existing waves, but we have not yet discussed how to detect the initial waves. Initial waves are ob­ tained around a given source point s E P . Let L be the layer containing s. The initial wave wa encloses a zero region containing this source point provided the region does not cover the complete layer L:

As an example, we interpret the first and second snap­ shots of figure 1. Snapshot 1 decomposes layer LI into two zero regions Au, A12 that are separated by a fron­ tier. The region Au contains the source points. Hence, Au is enclosed by the initial wave wo. wo cannot enclose A12 because otherwise the frontier between Au and AI2 is not explained. Therefore, Z1(wo) ={Au} .

In the second snapshot, layer £1 is divided into the zero regions A21 , An, A23, A24· The region A21 has fron­ tiers with An, A23, A24 . The old enclosed region of w0 is overlapping with A21 and An. Hence:

Region A22 is enclosed by w0 and touching the interface I 1 between the layers L1 and L3. Therefore, the incident e(wo, II) occurs in snapshot 2 and the waves r(e(wo, h)) and t( e( w0 , h)) exist in this snapshot. The reflected wave is in layer L 1 and its enclosed region is touching the interface h at the enclosed interface of the incident . Therefore, the reflected wave contains the region A22:

Layer L3 is divided into the zero regions B21 and B22 which are separated by a frontier. B21 is touching the interface h at the part that is touched by A22. Since the transmi tted wave is touching h at the enclosed interface of the incident we get:

In order to explain all observed frontiers, the waves can­ not enclose more regions as deduced above. Therefore , we get the following characterizations of waves:

This example shows that it is possible to interpret snap­ shots by a qualitative analysis of zero regions and their neighbourhood relationships. This method can be ex­ tended to other kinds of waves such as diffracted waves . Problems are encountered if 1. fronts end in a region without enclosing it and 2. two fronts of a different ori­ gin are linked without showing an indication where this link can be found . In order to treat these problems, we need addi t ional physical knowledge that cannot been ext racted from the images. In [Junker, 1994], a wave front has been characterized by a sequence of rays called polyray ( cf. figure 5) . Polyrays provide the additional knowledge, but are difficult to manage when traversing curved interfaces. A compromise could be the use of two auxiliary rays marking the left and right ends of wave fronts to meet the problems of the qualitative interpre­ tation method.

4 Related work

Recent work demonstrates the power of augmenting nu­ merical simulation with qualitative notions. [Forbus and Falkenhainer , 1990] define the notion of self-explanatory simulations where the simulator itself is able to explain its behavior :

a self-explanatory simulation integrates qualitative and numerical models to produce accurate predic­ tions and causal explanations of the behavior of continuous physical systems.

They illustrate this definition with the SIMGEN pro­ gram on physical systems simulated by ordinary differ­ ential equations.

Other examples of programs mixing quantitative sim­ ulation with qualitative notions can be found in the AI literature : the most famous being Q3 [Kuipers and Berleant, 1988] , POINCARE [Sacks, 1991], the Kineti­ cist's Workbench [Eisenberg, 1991], and others [Yip, 1987], [Zhao, 1991]. [Forbus, 1991] addressed an exten­ sion of qualitative reasoning to spatial information. In this work, Forbus advocates that, in order to be able to reason about spatially distributed system, one needs to mix two representations which he calls a metric dia­ gram (the quantitative part), and a place vocabulary (the symbolic part). The metl·ic diagram is used for calcula­ tion whereas the place vocabulary is used for describing the system's behavior at a more abstract level, and for guiding the numeric computations which take place on the metric diagram. Moreover, these two representations are intertwined so that there is a correspondence be­ tween the places identified by the place vocabulary, and the quantities manipulated in the analog representation. In a more recent paper [Forbus, 1994], he proposes six challenge problems for spatial reasoning, the fourth one being :

develop a system which can, given a sequence of weather maps for a region, provide a consistent qualitative explanation of the atmospheric behav­ ior cl uring that period ...

The problem we have addressed is very similar: Given a sequence of 2D snapshots of seismic amplitudes within the underground, our method provides a consistent qual­ itative explanation of the propagation of acoustic waves during that period . This has been achieved by effectively integrating several representations, namely a metric dia­ gram (i .e. fields) used for simulation and a place vocab­ ulary (i.e. objects) describing the geological structures. [Junker, 1994] additionally experimented with a physical representation based on rays.

Research in qualitative and model-based reasoning has focused since its beginning on systems that could be sim­ ulated by ordinary differential equations (ODEs). Nu­ merical simulators using differential equations can be di­ vided into two classes: Those using scalar variables and those using field variables. Scalar variables describe dif­ ferent quantitative properties of a system and are not distributed over a space. Good examples of this class of systems are simple physical devices, chemical pro­ cesses, chemical kinetics, global socio-economical mod­ els, or econometric models. In contrast to this, field variables are distributed over a space, often related to the real world in one, two or three dimensions. Within this category, we can distinguish between fields of scalar variables and fields of vector variables. Different simu­ lation techniques are used for approaching this kind of problems. Finite difference and finite elements are the conventional tools used by applied mathematicians for the simulation of field variables. Examples are fluid dy­ namics, geophysics or mechanics. Other approaches for field variables are naive physics and cellular automata, the basis for a number of ALife experiments such as Con­ way's game of life.

5 Conclusion

We developed a method for interpreting snapshot images produced by finite element simulators for seismic wave propagation. As a result, the regions in the images are linked with Hayes-like histories of waves:

1. In order to detect wave fronts in a snapshot, we characterized them by the zero regions in the back­ ground they are enclosing. The first snapshot con­ tains a single front enclosing the zero region that contains the source point.
2. Symbolic constraints are posed on the zero regions to track a given wave from one snapshot to the other and to detect new waves . We obtain new waves when wave histories intersect with the static histo­ ries of obstacles.

A first prototype of a snapshot interpreter which is called SISMONAUTE [Junker, 1994] has been implemented using the ILOG tools LELISP, AIDA, and SMECI. This expe­ rience enabled us to find the crucial concepts for char­ acterizing waves and for describing histories, as well as symbolic constraints, which enables the use of constraint programming tools to find globally consistent interpre­ tations.

As a future perspective, the interpretation method could be adapted to other kinds of numerical simulations (e.g. that of flame fronts in simulations of combustions).

Acknowledgements

This work was initiated by Alain Bamberger (IFP) and significantly supported by Laurence Nicoletis, Laurent Anne, Jean Brae (IFP) and Michel Gondran (EDF). Fur­ thermore, the authors would like to thank Claude Le Pape and the anonymous referees for their suggestions and comments.