WEBVTT
Kind: captions
Language: en

00:00:01.000 --> 00:00:13.640
[Silence]

00:00:13.640 --> 00:00:16.820
[background conversations]

00:00:16.820 --> 00:00:18.420
Good morning, everyone.

00:00:18.420 --> 00:00:21.660
Welcome to this week’s
Earthquake Science seminar.

00:00:22.180 --> 00:00:25.360
Before we start, I just want to
remind everyone, next week,

00:00:25.369 --> 00:00:27.310
Artie Rodgers will be here
from Lawrence Livermore

00:00:27.310 --> 00:00:29.740
talking about high-performance
computing of ground motion

00:00:29.740 --> 00:00:33.000
simulations – specifically about
the Hayward Fault scenario

00:00:33.000 --> 00:00:36.120
and some other smaller or
modest-magnitude events.

00:00:36.740 --> 00:00:41.820
Today – and, again, he just told me,
but I’m going to mispronounce this.

00:00:41.820 --> 00:00:48.100
Hejun Zhu from – who’s a assistant
professor of geophysics at UT-Dallas.

00:00:48.100 --> 00:00:51.239
Previously, he’s done a postdoc
at UT-Austin, where he held the

00:00:51.240 --> 00:00:53.940
Jackson School Distinguished
Postdoctoral Fellowship.

00:00:53.940 --> 00:01:00.260
And he, before that, completed his
Ph.D. with Jeroen Tromp at Princeton.

00:01:00.260 --> 00:01:04.750
His research has covered seismic imaging
of the Earth’s interior via tomography

00:01:04.750 --> 00:01:08.610
and migration techniques, which I think
we’re going to learn about today.

00:01:08.610 --> 00:01:12.250
Looking at uncertainty quantification
and geophysical inverse problems

00:01:12.250 --> 00:01:15.720
as well as some numerical simulations
of seismic wave propagation.

00:01:15.720 --> 00:01:17.820
Please join me in
welcoming our speaker.

00:01:17.820 --> 00:01:19.890
If you’d like to have
lunch with us afterwards,

00:01:19.890 --> 00:01:22.500
please find us after the seminar.
Or if you’d like to meet with

00:01:22.500 --> 00:01:27.300
Hejun before he takes off this afternoon
for SFO, please get a hold of me.

00:01:27.300 --> 00:01:30.600
- Okay, thanks, Alex.
Good morning, everyone.

00:01:30.600 --> 00:01:33.960
And I’m very happy to – here –
to share our perspectives

00:01:33.960 --> 00:01:39.210
towards elastic reverse time migration,
and also its applications on

00:01:39.210 --> 00:01:43.720
imaging Earth’s structures
as well as earthquake rupture.

00:01:43.720 --> 00:01:49.120
So this is a collaborative work
with my grad student, Jidong Yang.

00:01:49.970 --> 00:01:54.240
So first, please allow me to
introduce our research focus.

00:01:54.259 --> 00:01:58.100
We are computational
and structure seismologists.

00:01:58.100 --> 00:02:02.260
And we are interested in developing
advanced imaging technologies

00:02:02.260 --> 00:02:05.600
to investigate the physical properties
of the Earth’s materials.

00:02:05.600 --> 00:02:09.280
And we also want to
use those imaging technologies

00:02:09.280 --> 00:02:12.480
to characterize earthquake
rupture processes.

00:02:13.100 --> 00:02:16.680
So in my group, we are
working on these three directions.

00:02:16.680 --> 00:02:22.470
They are global seismology, exploration
seismology, and earthquake seismology.

00:02:22.470 --> 00:02:25.800
So to us, they are not
separated directions.

00:02:25.800 --> 00:02:28.850
They are actually internally
connected with each other.

00:02:28.850 --> 00:02:32.079
For instance, we all use
seismic wave equations

00:02:32.080 --> 00:02:35.260
to study seismic wave
propagation phenomena.

00:02:36.140 --> 00:02:41.519
So for instance, in global seismology,
we are working on a number of projects.

00:02:41.519 --> 00:02:45.549
On continent-scale tomography,
over the past several years,

00:02:45.549 --> 00:02:49.440
we have been constructing
3D mantle-scale tomographic models

00:02:49.440 --> 00:02:53.300
for the European and the
North American continents.

00:02:53.300 --> 00:02:57.260
And those results have
contributed to our current knowledge

00:02:57.260 --> 00:03:00.560
on the mantle dynamics
in those regions.

00:03:00.560 --> 00:03:04.220
And we have also been working on
crustal-scale tomography –

00:03:04.220 --> 00:03:08.400
for instance, right now,
we are focusing on two regions.

00:03:08.400 --> 00:03:14.700
One is north Texas and Oklahoma, and
the second one is the Delaware Basin.

00:03:14.700 --> 00:03:19.500
And hopefully, these 3D crustal
tomographic models can contribute to

00:03:19.500 --> 00:03:23.400
the current seismicity
studies in those regions.

00:03:23.400 --> 00:03:27.000
And we have also been working
on ambient noise tomography.

00:03:27.010 --> 00:03:31.660
And in terms of physical properties
of the Earth’s material, in particular,

00:03:31.660 --> 00:03:36.380
we are interested in azimuthal
anisotropy and shear attenuation,

00:03:36.380 --> 00:03:39.660
which can provide us
complementary information to

00:03:39.660 --> 00:03:44.440
the mantle dynamics in comparison
to seismic velocity distribution.

00:03:45.980 --> 00:03:50.160
In the exploration seismology side,
we are focusing on developing

00:03:50.160 --> 00:03:54.100
advanced imaging technologies.
So today we are working on

00:03:54.100 --> 00:03:58.600
three imaging technologies.
They are reverse time migration,

00:03:58.600 --> 00:04:01.870
full waveform inversion,
and Gaussian beam migration.

00:04:01.870 --> 00:04:06.450
This reverse time migration allows us
to better delineate subsurface

00:04:06.450 --> 00:04:11.130
distribution of reflectivities.
In its current stage, we are trying

00:04:11.130 --> 00:04:16.660
to promote elastic, anisotropic,
and attenuation RTM.

00:04:16.660 --> 00:04:21.380
Full waveform inversion is a imaging
technology, allows us to better constrain

00:04:21.380 --> 00:04:25.040
the subsurface velocity distributions.
So at the current stage, my group

00:04:25.040 --> 00:04:31.320
is tackling some technical issues in this
technology, including how to deal with

00:04:31.320 --> 00:04:36.590
a cycle-skipping problem, how to
mitigate multi-parameter tradeoffs,

00:04:36.590 --> 00:04:40.840
and also how to analyze
resolution and uncertainties.

00:04:40.840 --> 00:04:44.740
And this Gaussian beam migration
is one of ray-based method.

00:04:44.740 --> 00:04:47.090
In its current stage,
we are developing elastic

00:04:47.090 --> 00:04:50.380
and least square
Gaussian beam migration.

00:04:51.080 --> 00:04:55.420
In earthquake seismology side,
we are interested in locating earthquakes

00:04:55.420 --> 00:05:00.680
constraining moment tensor solutions
and imaging finite fault ruptures.

00:05:00.690 --> 00:05:06.070
So since we are structure seismologists,
our philosophy is Earth’s materials

00:05:06.070 --> 00:05:10.560
always involve lateral variation.
So we would like to – imaging – say,

00:05:10.560 --> 00:05:13.970
locate earthquakes,
constrain moment tensor solutions,

00:05:13.970 --> 00:05:18.400
and imaging this finite fault
rupture in fully 3D velocity models.

00:05:18.400 --> 00:05:23.460
So those are the three directions that we
are currently working on in my group.

00:05:23.460 --> 00:05:25.810
So let’s go to the
outline of today’s talk.

00:05:25.810 --> 00:05:32.160
First, I’m going to introduce to you the
basic idea of reverse time migration.

00:05:32.160 --> 00:05:36.700
We call it RTM – in case that you are not
familiar with exploration seismology.

00:05:36.700 --> 00:05:40.570
And then I will show our
developments on elastic

00:05:40.570 --> 00:05:44.940
reverse time migration and its
application on structure imaging.

00:05:44.940 --> 00:05:47.580
And after that,
we will show its application

00:05:47.580 --> 00:05:49.700
on earthquake
rupture imaging.

00:05:49.700 --> 00:05:54.930
In its current stage, we are trying to
use this technology on real data sets.

00:05:54.930 --> 00:05:57.930
So this is still
a work-in-progress.

00:05:57.930 --> 00:06:00.840
And finally, I will draw
some conclusions.

00:06:01.920 --> 00:06:06.220
Now, let’s first take
seismic migration 101.

00:06:06.220 --> 00:06:12.530
So for instance, this is a key
technology used in energy industry

00:06:12.530 --> 00:06:16.590
to find oil and gas reservoir.
For instance, if you have this very large

00:06:16.590 --> 00:06:23.170
vessel here which can tow a several-
kilometer streamer behind this vessel.

00:06:23.170 --> 00:06:27.420
And you excite this acoustic
energy by using those air guns.

00:06:27.420 --> 00:06:31.060
This will generate acoustic energy,
propagate through those water,

00:06:31.060 --> 00:06:35.810
and once it hits the impedance
contrast within the crust,

00:06:35.810 --> 00:06:40.170
those energy gets reflected and
recorded by the geophone attached to

00:06:40.170 --> 00:06:45.080
those streamers. And then we
can get this common-shot gather here.

00:06:45.080 --> 00:06:48.690
And in these gathers, you can
see those beautiful hyperbola.

00:06:48.690 --> 00:06:53.820
Those are reflection energy here.
So the seismic migration is trying to map

00:06:53.820 --> 00:06:58.030
those reflection energy from the time
domain to the depths domain so that

00:06:58.030 --> 00:07:02.920
we can better delineate the subsurface
distribution of the reflectivity.

00:07:03.889 --> 00:07:09.300
So over the past couple of decades,
there are lots of imaging technologies

00:07:09.310 --> 00:07:14.020
to work in the energy industry.
RTM is one of them.

00:07:14.020 --> 00:07:19.140
So this is probably – is the most
accurate depth imaging we have today,

00:07:19.140 --> 00:07:22.919
which enable us to map
the distribution of reflectivity

00:07:22.919 --> 00:07:25.530
in complicated
subsurface environments.

00:07:25.530 --> 00:07:30.230
So here I show you
three results from field data.

00:07:30.230 --> 00:07:33.040
So in this model,
we have this salt body,

00:07:33.040 --> 00:07:38.310
which has a very fast velocity.
And this probably is the

00:07:38.310 --> 00:07:43.030
common case in Gulf of Mexico.
So here I show you three panels.

00:07:43.030 --> 00:07:46.720
From left- to right-hand side,
we have migration results from

00:07:46.720 --> 00:07:51.760
Kirchhoff, one-way migration,
and the reverse time migration results.

00:07:51.770 --> 00:07:57.230
So I would like to draw your attention
to those events near the salt flank.

00:07:57.230 --> 00:08:00.990
So if we compare these three
technologies, you can see this RTM

00:08:00.990 --> 00:08:06.419
actually enable us to – delineates very
detailed structure near the salt flank.

00:08:06.419 --> 00:08:11.240
Because this is actually a critical zone
to accumulate oil and gas reservoir.

00:08:11.240 --> 00:08:15.240
So we needed to have
very detailed information here,

00:08:15.240 --> 00:08:21.260
where RTM enabled us to map very
detailed structure near the salt flank.

00:08:22.820 --> 00:08:26.060
So here I used – so these –
I borrowed this simple example to

00:08:26.060 --> 00:08:30.340
illustrate the basic framework, or basic
process, of reverse time migration.

00:08:30.340 --> 00:08:35.260
So there we have two steps.
And the first step records as a

00:08:35.260 --> 00:08:39.580
wavefield extrapolation, in which case,
we have to solve a wave equation

00:08:39.580 --> 00:08:42.120
in some complicated
velocity model.

00:08:42.120 --> 00:08:46.200
So typically, we have to
solve wave equation twice.

00:08:46.200 --> 00:08:51.310
In the first case, we put the source at
the source location, and then run the

00:08:51.310 --> 00:08:56.920
wavefield simulation forward to get
what we call it, as a source wavefield.

00:08:56.920 --> 00:09:02.020
In the second round, we inject recorded
data at surface and then run the wave

00:09:02.020 --> 00:09:06.620
equation backward to obtain what’s
recorded as a receiver wavefield.

00:09:06.620 --> 00:09:10.960
So those two steps are what we
call the wavefield extrapolation.

00:09:10.960 --> 00:09:13.990
And once you have
these two wavefields,

00:09:13.990 --> 00:09:18.980
we can apply some kind of imaging
condition to form migration imaging.

00:09:18.980 --> 00:09:24.860
For instance, this is the result if you
calculate the zero-lag cross-correlation

00:09:24.860 --> 00:09:29.850
between these two wavefields, we are
able to form this migration imaging.

00:09:29.850 --> 00:09:33.650
This is just for one shot.
Assuming that we have thousands

00:09:33.650 --> 00:09:37.130
or hundreds of shots,
we can stack these migration results

00:09:37.130 --> 00:09:41.360
together to get this result.
So you can see, in this result,

00:09:41.360 --> 00:09:45.820
we are able to delineate those
impedance contrast in the subsurface.

00:09:46.560 --> 00:09:49.740
So mathematically, we can write
down this imaging condition.

00:09:49.740 --> 00:09:55.380
For instance, we can obtain this
migration imaging as a function of space,

00:09:55.380 --> 00:09:58.930
which can be calculated by
using zero-lag cross-correlation

00:09:58.930 --> 00:10:03.480
between your source wavefield
and your receiver wavefield.

00:10:03.480 --> 00:10:07.870
So this is typically done in acoustic case,
therefore, we know that source

00:10:07.870 --> 00:10:12.340
wavefield and the receiver wavefield,
they are scalar fields.

00:10:13.040 --> 00:10:16.820
So how about if we want to
push it to elastic medium?

00:10:17.940 --> 00:10:21.370
The first, most straightforward,
one for extending these

00:10:21.370 --> 00:10:24.690
acoustic RTM to the elastic
medium is to use these

00:10:24.690 --> 00:10:27.900
component-wise cross-correlation
imaging condition.

00:10:27.900 --> 00:10:32.500
For instance, we can obtain this
vertical-vertical components migration

00:10:32.500 --> 00:10:36.360
imaging by using the zero-lag
cross-correlation between the vertical

00:10:36.360 --> 00:10:41.000
component source wavefield and the
vertical component’s receiver wavefield.

00:10:41.000 --> 00:10:43.780
So the same principle applies to
vertical-horizontal,

00:10:43.780 --> 00:10:47.360
horizontal-horizontal,
and horizontal-vertical.

00:10:47.360 --> 00:10:52.720
However, in this case, if your goal is
to delineate subsurface distribution

00:10:52.730 --> 00:10:58.070
of P-P or P-S reflectivity, obviously,
it is not a good choice to use

00:10:58.070 --> 00:11:02.751
these imaging conditions.
Because we know that P and S wave

00:11:02.751 --> 00:11:08.830
energies, they’re actually mixed in
vertical and horizontal components.

00:11:08.830 --> 00:11:13.620
So if you want to use elastic RTM,
the first step we would like to do

00:11:13.620 --> 00:11:18.780
is to decouple, or to separate,
compressional and shear wave energies.

00:11:18.780 --> 00:11:23.300
There are many ways to do that.
The simplest way is to apply this

00:11:23.300 --> 00:11:29.370
so-called divergence and curl separation.
For instance, I can obtain this P wave

00:11:29.370 --> 00:11:34.440
energy by applying this divergence
operator to your input wavefield.

00:11:34.440 --> 00:11:38.050
And then I can obtain
this S wave energy by applying

00:11:38.050 --> 00:11:41.529
this curl operator to
an input wavefield.

00:11:41.529 --> 00:11:48.980
So the thing to use these is that it's
very cheap and very easy to implement,

00:11:48.980 --> 00:11:52.240
but there are several issues.
For instance, in this case,

00:11:52.240 --> 00:11:55.460
the decoupled P wave is
a scalar field, and the

00:11:55.460 --> 00:11:58.190
decoupled S wave
is a vector field.

00:11:58.190 --> 00:12:01.810
But as we know, that in elastic
medium, both P and S wavefield,

00:12:01.810 --> 00:12:06.310
they should be vector fields.
And secondly, but using this

00:12:06.310 --> 00:12:10.520
spatial differential operator,
we actually introduce not only

00:12:10.520 --> 00:12:14.610
amplitude distortion, but also,
we introduce this phase shift.

00:12:14.610 --> 00:12:19.370
In other words, even the units of this
decoupled P waves is actually

00:12:19.370 --> 00:12:23.709
different from the input wavefield, u.
So we would like to design a better

00:12:23.709 --> 00:12:29.020
algorithm to better decouple the
vector P and S wave energies.

00:12:30.420 --> 00:12:35.899
And the second issue in elastic RTM is
that we have this polarity reversal issue.

00:12:35.899 --> 00:12:39.790
So if you remember your advanced
seismology class, you might remember

00:12:39.790 --> 00:12:45.540
that the P-S reflectivity coefficients are
actually – they’re not even functions.

00:12:45.540 --> 00:12:48.490
They’re actually odd function.
In which case, you are able to –

00:12:48.490 --> 00:12:52.330
you actually will find out that
the P-S reflectivity coefficient

00:12:52.330 --> 00:12:56.650
will change sign once you have
crossed the normal reflection point.

00:12:56.650 --> 00:13:01.100
So I will use this simple diagram to
illustrate that polarity reversal problem.

00:13:01.100 --> 00:13:05.480
So assuming that I have a source here,
and we have a tilted reflector,

00:13:05.480 --> 00:13:10.000
and the normal reflection point is here,
for instance, if I want to image in this

00:13:10.000 --> 00:13:14.320
A-prime location, by using the
zero-lag cross-correlation between

00:13:14.330 --> 00:13:18.520
the down-going P wave and
up-going S wave, or in this case,

00:13:18.520 --> 00:13:22.670
S wave actually pointing
towards negative x direction,

00:13:22.670 --> 00:13:26.670
then I’m going to assign
a negative sign at this location.

00:13:26.670 --> 00:13:31.020
While, in contrast, if you want to
image in this B-prime location,

00:13:31.020 --> 00:13:34.480
I can use the zero-lag cross-correlation
between this down-going P and

00:13:34.480 --> 00:13:38.300
up-going S, where in this case,
S actually pointing towards

00:13:38.300 --> 00:13:41.779
positive x direction.
Therefore, I actually assign

00:13:41.780 --> 00:13:43.640
a plus sign
at this location.

00:13:43.640 --> 00:13:45.140
In other words,
when you cross the

00:13:45.150 --> 00:13:49.499
normal reflection point, you have
this polarity reversal issue.

00:13:49.499 --> 00:13:53.150
And you will see later, in a numerical
example, you will see this is actually a

00:13:53.150 --> 00:13:59.600
very serious issue when you want to get
very accurate P-S migration images.

00:13:59.600 --> 00:14:04.430
So we would like to tackle these two
issues by using better vector P and

00:14:04.430 --> 00:14:09.040
S wave decomposition algorithms.
So what we did is, we go back to

00:14:09.040 --> 00:14:13.200
Aki and Richards' classical
Quantitative Seismology.

00:14:13.200 --> 00:14:16.550
They actually provided us
a recipe to better decompose

00:14:16.550 --> 00:14:18.930
P and S wave energies.

00:14:18.930 --> 00:14:23.320
The tool that we are going to use is
called the Helmholtz decomposition.

00:14:23.320 --> 00:14:26.830
For instance, I know the
wavefield, u, which include

00:14:26.830 --> 00:14:34.100
both vector P and S energy here.
And here I introduce two potentials.

00:14:34.100 --> 00:14:38.740
One is a scalar potential, phi.
And one is a vector potential, psi.

00:14:38.740 --> 00:14:42.380
And the scalar potential, phi,
is related to P wave energy.

00:14:42.390 --> 00:14:47.480
And the vector potential, psi,
is related to S wave energy.

00:14:47.480 --> 00:14:51.290
So the question is, let’s –
if I give you those two potentials,

00:14:51.290 --> 00:14:57.029
I can easily calculate P and S wave
energy and also the input wavefield.

00:14:57.029 --> 00:14:59.940
But right now, the question
is the other way around, is that,

00:14:59.940 --> 00:15:03.089
if I give you this wavefield, u,
how can you calculate

00:15:03.089 --> 00:15:06.200
these scalar and
vector potentials?

00:15:07.080 --> 00:15:10.839
And Aki and Richards actually
provided us the way to do that,

00:15:10.839 --> 00:15:15.260
is let’s – we actually need to
introduce this additional wavefield, w,

00:15:15.260 --> 00:15:19.310
which has to satisfy this
vector Poisson equation.

00:15:19.310 --> 00:15:23.080
And if you remember this identity
relationship with a Laplacian

00:15:23.080 --> 00:15:27.500
operator, by comparing these three
equations, you might immediately

00:15:27.510 --> 00:15:31.850
recognize that this scalar potential
can be calculated by using

00:15:31.850 --> 00:15:34.470
the divergence of this
additional wavefield, w.

00:15:34.470 --> 00:15:37.390
And this vector potential
can be calculated by using

00:15:37.390 --> 00:15:41.180
the [inaudible] curl of this
additional wavefield, w.

00:15:41.180 --> 00:15:44.899
And then once I solve this w,
P and S wave energies

00:15:44.900 --> 00:15:47.980
can be decoupled by
using these two equations.

00:15:47.980 --> 00:15:52.480
So the key question is, how do
you solve this w wavefield?

00:15:53.920 --> 00:15:58.680
In the classical Helmholtz
decomposition, if you look at this

00:15:58.680 --> 00:16:02.730
vector Poisson’s equation, we do
have a Green’s function for this case.

00:16:02.730 --> 00:16:07.300
Because this Poisson equation
has constant coefficients,

00:16:07.300 --> 00:16:09.950
so we have a Green’s
function, look like this.

00:16:09.950 --> 00:16:12.750
And once you know the Green’s
function by using the representation

00:16:12.750 --> 00:16:16.980
theorem, we can directly write down
an analytical solution to decouple

00:16:16.980 --> 00:16:19.380
P and S wave energy,
shown here.

00:16:19.380 --> 00:16:22.579
All right, the equations,
they are very simple.

00:16:22.579 --> 00:16:25.860
But they’re actually useless
because they’re actually

00:16:25.860 --> 00:16:31.840
very computationally involved.
Because in order to calculate those

00:16:31.840 --> 00:16:36.860
3D volume integrals, the computational
complexity is the order of N-squared.

00:16:36.860 --> 00:16:41.680
So if we are dealing with 3D problems,
this N can go to 10 to the 9th.

00:16:41.680 --> 00:16:45.100
Which means, in order to evaluate
those 3D volume integrals,

00:16:45.100 --> 00:16:48.940
the computational complexity
can go to 10 to the 18.

00:16:49.500 --> 00:16:55.040
So we have to find an efficient algorithm
to solve this vector Poisson’s equation.

00:16:55.040 --> 00:16:58.910
And the good news is that in
numerical analysis, we do have

00:16:58.910 --> 00:17:03.470
a number of efficient algorithms to
solve this vector Poisson’s equation.

00:17:03.470 --> 00:17:06.920
The first Poisson solver
is one of these algorithm.

00:17:06.920 --> 00:17:10.900
So there are three steps
in the fast Poisson solver.

00:17:10.910 --> 00:17:14.390
In the first step, we calculate
the Fourier spectra of your

00:17:14.390 --> 00:17:19.130
input wavefield, U,
by using the forward sine transform.

00:17:19.130 --> 00:17:21.989
And in the second step,
we calculate the Fourier spectra

00:17:21.989 --> 00:17:26.610
of this W wavefield by using
this element-wise division.

00:17:26.610 --> 00:17:31.850
And once you know the Fourier spectra
of this W, I can group it and get it

00:17:31.850 --> 00:17:35.240
back to the physical domain by
using inverse sine transform.

00:17:35.240 --> 00:17:39.289
So the reason we can use those
sine function and this specially

00:17:39.289 --> 00:17:44.279
chosen lambda coefficients is because
those sine functions and the lambda,

00:17:44.279 --> 00:17:48.049
they are actually eigenvalue
and eigenvectors of this

00:17:48.049 --> 00:17:50.220
discretized
Laplacian operator.

00:17:50.220 --> 00:17:54.580
So I can use this algorithm to
solve this Poisson’s equation.

00:17:54.590 --> 00:18:00.519
And once you solve this W wavefield,
I can calculate this vector P and S wave.

00:18:00.519 --> 00:18:05.169
So the advantages that we are able to
obtain are vector P and S wavefields.

00:18:05.169 --> 00:18:08.680
And you have correct
amplitude as well as the phase.

00:18:08.680 --> 00:18:11.769
Now, let’s count the
computational complexity.

00:18:11.769 --> 00:18:16.859
For the step 1 and the step 3, we can use
the fast Fourier transform to solve this

00:18:16.860 --> 00:18:20.720
sine transform, so the computational
complexity is the order of N-log-N.

00:18:20.720 --> 00:18:24.460
And step 1 – step 2,
just element-wise division.

00:18:24.460 --> 00:18:28.440
So if you count the total computational
complexity, it’s just N-log-N,

00:18:28.450 --> 00:18:32.610
which is much more efficient than
the previous analytical solution.

00:18:32.610 --> 00:18:34.289
So in this study,
we actually utilized this

00:18:34.289 --> 00:18:39.540
fast Poisson solver to better
decouple P and S wave energies.

00:18:40.560 --> 00:18:45.080
So now let’s test this algorithm for a
simple homogeneous velocity model.

00:18:45.080 --> 00:18:49.619
So for instance, this is one snapshot
of vertical components and horizontal

00:18:49.619 --> 00:18:56.249
components displacements, which
involve both P wave and S wave here.

00:18:56.249 --> 00:18:59.780
If I just apply the divergence
and curl separation,

00:18:59.780 --> 00:19:02.809
those are the potentials
that I’m able to obtain.

00:19:02.809 --> 00:19:06.320
And if you look carefully, for instance,
if you look at this scale bar,

00:19:06.320 --> 00:19:10.030
they’re actually totally
different from our input wavefield.

00:19:10.030 --> 00:19:13.809
And also, if you look for those phases,
because you are using those

00:19:13.809 --> 00:19:18.490
differential operator, as you actually 
introduce those phase shift,

00:19:18.490 --> 00:19:22.600
if you compare this potential
and this P wave energy here.

00:19:23.420 --> 00:19:27.900
So if you apply our Helmholtz
decomposition to that wavefield,

00:19:27.900 --> 00:19:33.040
those are the decoupled P waves on
vertical and horizontal components.

00:19:33.040 --> 00:19:37.740
There’s also the decoupled S wave
on vertical and horizontal components.

00:19:37.740 --> 00:19:43.379
So the beauty of this algorithm is that, if
I sum over this P wave and S wave along

00:19:43.379 --> 00:19:49.490
horizontal direction, I’m able to perfectly
recover the input horizontal wavefields.

00:19:49.490 --> 00:19:54.779
So this is just for simple velocity model.
Now, let’s test this algorithm for

00:19:54.779 --> 00:19:59.720
more complicated velocity model. So
this is a classical bemchmark model –

00:19:59.720 --> 00:20:04.330
3D SEG/EAGE salt model,
in which case we have this very fast

00:20:04.330 --> 00:20:09.190
velocity heterogeneity of salt
body in the middle of your model.

00:20:09.190 --> 00:20:15.279
And if I solve 3D elastic wavefields,
this is one snapshot we have along

00:20:15.279 --> 00:20:19.809
X, Y, and Z directions,
in which case, we have both

00:20:19.809 --> 00:20:24.640
compressional and shear
energy in those snapshots.

00:20:24.640 --> 00:20:28.440
And in the second row,
I show you decoupled P wave

00:20:28.440 --> 00:20:32.009
along X, Y,
and Z directions.

00:20:32.009 --> 00:20:37.429
And in this case, we only have
compressional energy in those snapshots.

00:20:37.429 --> 00:20:43.480
And the third row shows decoupled
S waves around X, Y, and Z directions.

00:20:43.480 --> 00:20:49.559
So the beauty is that, if I sum always a
P wave and S wave along the Y direction,

00:20:49.560 --> 00:20:55.380
I’m able to perfectly recover the
inputs wavefield around the Y axis.

00:20:56.660 --> 00:21:01.700
So by using this algorithm, I can
decouple vector P and S wave energies.

00:21:01.700 --> 00:21:05.200
In this case, I can introduce this
vector imaging condition to

00:21:05.200 --> 00:21:09.029
better form migration results
when we perform elastic RTM.

00:21:09.029 --> 00:21:13.139
For instance, at this time,
I can form a P-P migration imaging

00:21:13.139 --> 00:21:17.999
by using the dot products
between the vector P wave in the

00:21:17.999 --> 00:21:21.200
source side and the vector
P wave in the receiver side.

00:21:21.200 --> 00:21:28.200
So the same principle applies to S-S,
P-S, and S-P migration imagings.

00:21:28.200 --> 00:21:32.720
Now let’s first test this for this relatively
simple velocity model, in which case,

00:21:32.720 --> 00:21:35.640
you have these
horizontal layers.

00:21:35.640 --> 00:21:40.489
And this is the migration imaging if I use
vertical-vertical components migration.

00:21:40.489 --> 00:21:44.960
So the imaging quality is fine.
But, as we mentioned, the P and S

00:21:44.960 --> 00:21:48.720
waves, they are actually coupled in
this vertical and horizontal component,

00:21:48.720 --> 00:21:52.240
therefore, we do not have
a specific meaning – physical meaning

00:21:52.240 --> 00:21:54.800
for this migration imaging.
And if you just use

00:21:54.800 --> 00:21:59.460
divergence-to-curl separation,
this is the P-P results we have.

00:21:59.460 --> 00:22:01.070
And if you use the
Helmholtz decomposition

00:22:01.070 --> 00:22:04.710
and the vector imaging condition,
this is the P-P results we have.

00:22:04.710 --> 00:22:09.440
Well, if you look at the P-P,
you might say the image [inaudible]

00:22:09.440 --> 00:22:12.149
by using this
divergence-to-curl is fine.

00:22:12.149 --> 00:22:14.940
You don’t see any
significant improvement.

00:22:14.940 --> 00:22:18.760
While, if you look at the P-S imaging,
for instance, if you use

00:22:18.760 --> 00:22:22.380
divergence-to-curl, and if you do
not apply any sophisticated

00:22:22.380 --> 00:22:27.340
angle-dependent correction terms,
this is the migration result you have.

00:22:27.340 --> 00:22:30.059
And I would like to draw
your attention for this reflector.

00:22:30.059 --> 00:22:34.580
For instance, when you cross
these points, the color actually

00:22:34.580 --> 00:22:40.320
change from black, white, black
to white, black, white.

00:22:40.320 --> 00:22:44.639
So this is exactly the polarity
reversal issue I show you before.

00:22:44.639 --> 00:22:48.059
And if you apply the Helmholtz
decomposition and vector imaging

00:22:48.059 --> 00:22:52.289
condition, we’re actually able to
avoid this polarity reversal issue.

00:22:52.289 --> 00:22:56.429
In this case, we are able to
get much better reflector here.

00:22:56.429 --> 00:22:59.980
All right, so this is for this
simple velocity model.

00:22:59.980 --> 00:23:03.580
And later, my grad students tried
this for a Marmousi model,

00:23:03.580 --> 00:23:08.460
which is a complicated velocity model.
This is the P-S migration imaging if

00:23:08.460 --> 00:23:11.559
you just use divergence-to-curl
separation and if you do not

00:23:11.559 --> 00:23:17.029
apply any correction term. So we can see
this imaging quality is actually very bad.

00:23:17.029 --> 00:23:21.260
It’s very difficult for you
to interpret some key features.

00:23:21.260 --> 00:23:24.390
While, if you use the Helmholtz
decomposition and vector imaging

00:23:24.390 --> 00:23:27.739
condition, this is the
P-S migration results we have.

00:23:27.739 --> 00:23:31.059
We are able to
get much better P-S results

00:23:31.060 --> 00:23:34.680
in comparison to
the previous studies.

00:23:36.149 --> 00:23:38.480
So this technology we –
at current stage,

00:23:38.480 --> 00:23:42.020
we just test this for a
reservoir-scale imaging.

00:23:42.020 --> 00:23:47.179
And in the next step, we’d also like
to use it for mantle-scale studies.

00:23:47.179 --> 00:23:52.049
For instance, we would like to use this
technology to image discontinuities,

00:23:52.049 --> 00:23:57.340
say to cross the Moho or, like,
410 and 660 discontinuities,

00:23:57.340 --> 00:24:01.740
and we would also like to
image scatters within the mantle.

00:24:01.759 --> 00:24:05.190
The requirements to apply this
technology is that we need to

00:24:05.190 --> 00:24:08.799
have a really decent velocity
model in the subsurface.

00:24:08.799 --> 00:24:13.340
The good news is that, in my group,
over the past couple of years,

00:24:13.340 --> 00:24:18.429
we are trying to construct very decent
tomographic models for several regions.

00:24:18.429 --> 00:24:22.779
For instance, this is the current
tomographic model we have for the

00:24:22.779 --> 00:24:26.769
European continent by using
180 earthquakes and

00:24:26.769 --> 00:24:31.629
about 700 seismic stations
by using the advanced full waveform

00:24:31.629 --> 00:24:35.599
inversion technology to
build this velocity model.

00:24:35.599 --> 00:24:39.799
We also have a new tomographic
model for North American continent

00:24:39.799 --> 00:24:43.309
by using U.S. array data.
So hopefully, we would like to use

00:24:43.309 --> 00:24:49.119
this imaging technology to better
delineate, say, 410 and 660 or to find

00:24:49.120 --> 00:24:54.500
scatters within the mantle by using
those velocity tomographic models.

00:24:55.629 --> 00:25:00.640
So next, I would like to show, how do
we apply this elastic RTM to better

00:25:00.640 --> 00:25:05.120
image earthquakes? Especially
the earthquake rupture processes.

00:25:05.720 --> 00:25:11.160
So let’s go back to this imaging principle
proposed by Claerbout in 1971,

00:25:11.160 --> 00:25:16.380
where he said those reflectors
exist at the point where your

00:25:16.400 --> 00:25:20.400
down-going source wavefield
and up-going receiver wavefield,

00:25:20.400 --> 00:25:22.520
they are correlated
with each other.

00:25:22.520 --> 00:25:26.500
So mathematically, we can
use this previously introduced

00:25:26.500 --> 00:25:31.220
dot product imaging condition to
image those reflectors, for instance,

00:25:31.229 --> 00:25:32.989
by using the zero-lag
cross-correlation between

00:25:32.989 --> 00:25:37.409
P and S to imaging
the location of this reflector.

00:25:37.409 --> 00:25:42.360
However, if you want to use this idea –
this concept to image earthquakes,

00:25:42.360 --> 00:25:46.460
we actually get in trouble because,
in this case, we don’t have this

00:25:46.479 --> 00:25:50.389
down-going source wavefield.
We only have this receiver wavefield.

00:25:50.389 --> 00:25:52.820
So once you have
this receiver wavefield,

00:25:52.820 --> 00:25:56.340
how are you going to map,
say, the earthquake here?

00:25:56.340 --> 00:26:02.560
So we’re trying to do that by leveraging
the P and S wave separation in elastic

00:26:02.560 --> 00:26:06.599
medium to better locate earthquakes.
For instance, if you look at –

00:26:06.599 --> 00:26:10.909
for this homogeneous velocity model,
if you look at the receiver side,

00:26:10.909 --> 00:26:16.470
we have this first-arrival P and
a second arrival S wave energy.

00:26:16.470 --> 00:26:21.249
So if I am able to separate P and S
wave energy, and then compute

00:26:21.249 --> 00:26:26.700
zero-lag cross-correlation between these
two traces, I actually get nothing here.

00:26:27.580 --> 00:26:29.980
The same principle
applies to this location.

00:26:30.740 --> 00:26:33.760
However, when you’re near –
or, when you’re in the vicinity

00:26:33.769 --> 00:26:37.470
of the source location,
and if you decouple P and S wave,

00:26:37.470 --> 00:26:39.519
and then compute the
zero-lag cross-correlation

00:26:39.519 --> 00:26:44.250
between these two traces,
you actually will get some energy here.

00:26:44.250 --> 00:26:48.500
But if you move away from the
source location, you still get nothing.

00:26:48.500 --> 00:26:53.650
So that actually – inspire us to design
a imaging condition look like this one,

00:26:53.650 --> 00:26:57.240
where, in this case, instead of using
the dot product between the source

00:26:57.240 --> 00:26:59.999
and receiver wavefield,
we use the dot product between

00:26:59.999 --> 00:27:07.119
the P wave in the receiver side and
S wave on – also on the receiver side

00:27:07.119 --> 00:27:12.470
to form this migration imaging,
which allows us to locate earthquakes.

00:27:12.470 --> 00:27:17.580
So instead of using this summation
over the entire time duration,

00:27:17.580 --> 00:27:20.789
we also propose to use this
moving window summation.

00:27:20.789 --> 00:27:25.820
So, in this case, we are actually able
to not only locate a single earthquake,

00:27:25.820 --> 00:27:31.440
but to capture the entire rupture
processes, or fracture processes.

00:27:31.440 --> 00:27:36.740
So first, we would like to test our
algorithm by using a [inaudible] model.

00:27:36.740 --> 00:27:42.380
So this, again, is the Marmousi model.
So we have one single earthquake here.

00:27:42.380 --> 00:27:44.909
And assuming that I have
a whole bunch of stations

00:27:44.909 --> 00:27:51.440
at Earth surface, those are the recorded
X and Z component displacements.

00:27:52.100 --> 00:27:56.280
And if I use the imaging condition
proposed in the previous slide, those are

00:27:56.289 --> 00:28:00.499
the imaging that I’m able to obtain.
So here I would like to compare

00:28:00.500 --> 00:28:04.980
what’s happening if you use a
inaccurate P and S wave separation.

00:28:04.980 --> 00:28:08.980
So for instance, this is the imaging
we have if you use divergence-to-curl

00:28:08.980 --> 00:28:14.499
separation. And this is the case if
we use vector imaging condition.

00:28:14.499 --> 00:28:18.739
So you can see, from this view,
you don’t see any large differences.

00:28:18.740 --> 00:28:23.620
But if I zoom in, if you do not use
very accurate P and S wave separation,

00:28:23.620 --> 00:28:27.120
in this case, you actually locate
two earthquakes instead of one.

00:28:27.120 --> 00:28:31.279
While if you’re using our proposed
vector imaging condition,

00:28:31.279 --> 00:28:34.240
we do – able to narrow
down to a single event.

00:28:34.240 --> 00:28:38.640
And those spreading actually
tells you the resolution we have.

00:28:38.640 --> 00:28:43.309
So this is the case where you have very
clean data. You don’t have any noise.

00:28:43.309 --> 00:28:46.659
What happen if I applied
a Gaussian noise to my data?

00:28:46.659 --> 00:28:50.239
So here we test the three scenarios.
The first is a signal-to-noise

00:28:50.240 --> 00:28:56.940
ratio equals to 10, 5, and 1.
Before we back-project those

00:28:56.940 --> 00:29:01.760
recorded data to the subsurface,
those are the migration results we have.

00:29:01.760 --> 00:29:06.400
So for instance, for signal-to-noise
ratio equals to 10 and 5, we’re able to –

00:29:06.400 --> 00:29:09.860
still able to get those
earthquake locations.

00:29:09.860 --> 00:29:13.759
But, see if we have signal-to-noise
ratio equals to 1, we’re still able to

00:29:13.759 --> 00:29:17.259
get this location, but you suffer
from those low wave number

00:29:17.259 --> 00:29:21.240
artifacts at shallow depths.
So at current state, we are still trying

00:29:21.240 --> 00:29:27.299
to understand the reason for
those low wave number artifacts.

00:29:27.299 --> 00:29:31.750
So this is the case where
we have very dense coverage.

00:29:31.750 --> 00:29:35.600
So – oh, by the way, in this case,
you can see – if you look at this data set,

00:29:35.600 --> 00:29:40.900
there is no way that you
can pick any effective arrivals.

00:29:40.900 --> 00:29:45.320
So in our algorithm, we actually
do not need to pick any arrivals.

00:29:45.320 --> 00:29:48.500
So what happen if I don’t have
very good coverage, in which case –

00:29:48.500 --> 00:29:53.679
for instance, here, eliminates 20%
data and only preserve 80% data?

00:29:53.679 --> 00:29:57.759
Now, let’s back-project
those data sets to the subsurface.

00:29:57.759 --> 00:30:00.009
So this is, again, if you use
vector imaging condition,

00:30:00.009 --> 00:30:05.749
we’re still able to locate this earthquake.
And for signal-to-noise ratio equals to 1,

00:30:05.749 --> 00:30:09.809
we still suffer from those low wave
number artifacts at shallow depths,

00:30:09.809 --> 00:30:14.440
but you’re still being able to locate
these events at greater depths.

00:30:16.090 --> 00:30:18.860
So this is the case where we
can use the imaging condition

00:30:18.879 --> 00:30:22.049
to locate a single event.
So, as I promised, that we would

00:30:22.049 --> 00:30:26.600
also like to use it to image the
entire earthquake rupture processes.

00:30:26.600 --> 00:30:32.700
So we did this simple numerical test,
is that we used this 3D SEG/EAGE

00:30:32.700 --> 00:30:37.730
overthrust model, in which case, we have
this very strong – lot of heterogeneity.

00:30:37.730 --> 00:30:41.380
And then we put a vertical fault
plane in the middle of the model.

00:30:41.380 --> 00:30:46.049
So these two panels,
they are the vertical cross-section.

00:30:46.049 --> 00:30:50.300
One is along the fault plane, and one
is perpendicular to the fault plane.

00:30:50.300 --> 00:30:53.259
And this is the map of
the [inaudible] fault plane.

00:30:53.259 --> 00:30:58.109
And we discretized this fault plane
as whole bunch of grid points.

00:30:58.109 --> 00:31:02.169
And each grid points were
assigned different slip magnitude

00:31:02.169 --> 00:31:06.129
and also different initial time.
For instance, here, the color bar

00:31:06.129 --> 00:31:09.340
means the slip magnitude.
Therefore, you can see here

00:31:09.340 --> 00:31:15.580
we have a major slip patch occurred
at time less than 0.5 second.

00:31:15.580 --> 00:31:19.690
And a secondary slip patch
occurred at 1.5 second.

00:31:19.690 --> 00:31:24.989
And a third one was much weaker
amplitude occurred at 2.5 second.

00:31:24.989 --> 00:31:30.960
And those red dashed lines showed
initial time of each grid point.

00:31:30.960 --> 00:31:34.349
So our goal is actually to
image the entire spatial

00:31:34.349 --> 00:31:38.300
and temporal evolution
of these rupture processes.

00:31:39.240 --> 00:31:41.960
So let’s see how our migration
imaging look like.

00:31:41.969 --> 00:31:47.269
So this is typically how we present
a 3D volume on a 2D plane.

00:31:47.269 --> 00:31:51.950
So for instance, this is a vertical
plane along the fault’s plane.

00:31:51.950 --> 00:31:57.840
And this side view is the vertical plane
perpendicular to the – to the fault plane.

00:31:57.840 --> 00:32:02.259
And top one is actually a map
view across the fault plane.

00:32:02.259 --> 00:32:06.450
So on the left-bottom corner,
I show you the [inaudible] inputs –

00:32:06.450 --> 00:32:08.320
the slip distribution.

00:32:08.320 --> 00:32:11.690
So this is the migration imaging
we have as time equals to zero.

00:32:11.690 --> 00:32:17.300
We are able to capture this major slip
patch in the middle of the fault plane.

00:32:17.300 --> 00:32:22.639
And, as time evolve, we can see this –
we can image these rupture processes.

00:32:22.639 --> 00:32:25.889
For instance, even we are
able to imaging this

00:32:25.889 --> 00:32:30.640
secondary slip patch
occurred at 1.5 second here.

00:32:31.600 --> 00:32:36.580
And even this third one,
occurred at 2.5 second.

00:32:38.280 --> 00:32:42.039
So you can see this technology
actually allows us to locate earthquakes

00:32:42.039 --> 00:32:45.409
even to – in the future,
probably to monitor microseismicity

00:32:45.409 --> 00:32:48.899
and to imaging the fracture
or rupture propagation.

00:32:48.899 --> 00:32:52.769
So the advantage of using this
technology is that we do not need to

00:32:52.769 --> 00:32:55.900
invoke any 1D
velocity distributions.

00:32:55.900 --> 00:33:00.620
We can imaging this earthquake
by using fully 3D velocity models.

00:33:00.620 --> 00:33:05.460
And the fault actually does not
need to be vertical – a single plane.

00:33:05.460 --> 00:33:11.640
It can have a complex geometry, and
we do not need to apply phase picking.

00:33:13.180 --> 00:33:17.659
So if we compare the slip
distribution we have for the inputs,

00:33:17.659 --> 00:33:20.789
and the migrated results –
the cumulative slip distribution we have,

00:33:20.789 --> 00:33:24.529
and we’ve compared these two.
You can see we are able to image in

00:33:24.529 --> 00:33:29.830
this major slip patch, the secondary
slip patch, and the third one here.

00:33:29.830 --> 00:33:32.759
The third one, the magnitude
is relatively weak

00:33:32.760 --> 00:33:37.120
because we lose
resolution at greater depths.

00:33:38.380 --> 00:33:41.830
So we’d like to – in the next step,
we’d like to apply this technology

00:33:41.830 --> 00:33:45.440
to real data sets.
So this is still work in progress.

00:33:46.289 --> 00:33:50.080
So in order to doing so –
in order to do so, we actually have to

00:33:50.080 --> 00:33:53.660
have very decent 3D
crustal velocity models.

00:33:53.660 --> 00:33:57.960
So we would like to try our
technologies for two sites.

00:33:57.979 --> 00:34:02.159
The first one is for southern California,
in which case, we have very decent

00:34:02.160 --> 00:34:06.240
3D crustal velocity models,
and we have very good data coverage.

00:34:06.240 --> 00:34:11.000
And the second site we would
like to try is for central Oklahoma.

00:34:11.800 --> 00:34:14.860
So the reason we’d like to focus
on central Oklahoma is,

00:34:14.869 --> 00:34:19.159
over the past couple of years,
we observed significant increase

00:34:19.159 --> 00:34:23.000
in terms of seismicity in those area.
For instance, this diagram show the

00:34:23.000 --> 00:34:27.410
epicentral location in Texas
and Oklahoma, and you compare

00:34:27.410 --> 00:34:30.560
the number of earthquakes
between Oklahoma and California,

00:34:30.560 --> 00:34:35.060
you do observe significant
increase of seismicity.

00:34:35.060 --> 00:34:39.050
And we do – we did have
nine earthquakes with magnitude

00:34:39.050 --> 00:34:44.040
greater than 4.5 from 2011 to 2016.
But we would like to apply our

00:34:44.040 --> 00:34:48.450
technology to better investigate those
earthquakes and also, for instance,

00:34:48.450 --> 00:34:53.160
to better imaging the rupture
processes for those large earthquakes.

00:34:53.960 --> 00:34:57.420
So, as I said, if we want to use
this technology, first of all,

00:34:57.430 --> 00:35:00.380
we have to first image
the velocity crustal structure.

00:35:00.380 --> 00:35:03.140
So this current stage,
we are trying to do this.

00:35:03.140 --> 00:35:07.070
So in order to do that, we’re actually
trying to use earthquake tomography

00:35:07.070 --> 00:35:11.480
to better constrain the crustal
velocity model in central Oklahoma.

00:35:11.480 --> 00:35:16.340
So, in this project, we actually collect
more than 200 earthquakes and more

00:35:16.340 --> 00:35:21.020
than 200 stations in central Oklahoma.
And those earthquakes with

00:35:21.020 --> 00:35:26.960
magnitude range from 3 to 5, and
we stop shallower than 10 kilometer.

00:35:28.470 --> 00:35:30.880
So the knowledge that we are
going to use to image the

00:35:30.880 --> 00:35:34.450
crustal velocity distribution
is called full waveform inversion.

00:35:34.450 --> 00:35:37.090
So we’d like to leverage
the entire wavefield

00:35:37.090 --> 00:35:40.000
to better constrain
3D velocity structures.

00:35:40.000 --> 00:35:45.790
So for instance, we test 1D simulation.
We discretize – use a spectral-element

00:35:45.790 --> 00:35:53.040
mesh to discretize a 200-kilometer-by-
200-kilometer-by-50-kilometer box.

00:35:53.040 --> 00:35:58.680
And at first round, we just use a
1D velocity model from St. Louis

00:35:58.680 --> 00:36:04.170
moment tensor solution catalog. And
this is 1D velocity distribution we have.

00:36:04.170 --> 00:36:07.860
And we can check the comparisons
between data and the prediction

00:36:07.860 --> 00:36:12.170
from that 1D velocity model
for a single earthquake here.

00:36:12.170 --> 00:36:15.560
So on the right-hand side, from left-hand
side to right-hand side, we have vertical

00:36:15.560 --> 00:36:19.760
components, radial components, and the
transverse component seismograms.

00:36:19.760 --> 00:36:22.780
And then [inaudible]
as a function of distance.

00:36:22.780 --> 00:36:26.800
And the red seismogram,
the simulated seismogram

00:36:26.800 --> 00:36:32.250
for 1D velocity model, and the black
seismograms, they are real observations.

00:36:32.250 --> 00:36:35.760
So if you look at this
comparison, first of all,

00:36:35.760 --> 00:36:40.210
we do have a very good
fitting for P and S wave arrivals.

00:36:40.210 --> 00:36:44.510
But if you look at those [inaudible]
surface waves, obviously this 1D

00:36:44.510 --> 00:36:49.350
velocity model cannot capture these
complicated surface waves arrivals.

00:36:49.350 --> 00:36:53.140
And if I measure the travel time
differences between the observed

00:36:53.140 --> 00:36:56.790
and the predicted seismograms,
you might see that we do have

00:36:56.790 --> 00:37:00.340
a systematic bias in this case.
For instance, especially for radial

00:37:00.340 --> 00:37:04.140
and the transverse components,
the mean travel time differences

00:37:04.140 --> 00:37:08.590
is about a negative 5 second,
which suggests that 1D starting model

00:37:08.590 --> 00:37:11.590
is a little bit slow,
and we need to speed it up.

00:37:11.590 --> 00:37:15.200
So our goal is to actually
correct the systematic bias

00:37:15.200 --> 00:37:18.660
and reduce this
standard deviation.

00:37:18.660 --> 00:37:20.740
And, as I mentioned,
the technology that we are

00:37:20.750 --> 00:37:23.370
going to use is called
full waveform inversion.

00:37:23.370 --> 00:37:28.030
So I would like to use a couple of slides
to introduce why we would like to use

00:37:28.030 --> 00:37:32.640
that full waveform inversion technology
to constrain the velocity distributions.

00:37:32.640 --> 00:37:35.460
So this is how we do a
classical travel time tomography.

00:37:35.460 --> 00:37:40.750
For instance, if you have data sets like
this, first, you pick the arrival time,

00:37:40.750 --> 00:37:44.570
presumably by using first arrival.
And it puts this observed arrival

00:37:44.570 --> 00:37:48.590
time and this offset
in the travel time diagram.

00:37:48.590 --> 00:37:53.670
The secondary ones, you assume that
you know some a priori velocity model.

00:37:53.670 --> 00:37:58.040
Then we can trace seismic array
through this smooth velocity model

00:37:58.040 --> 00:38:02.940
and predict the travel time
as shown by this solid white line.

00:38:02.940 --> 00:38:06.410
And our goal is to minimize the
differences between these observed

00:38:06.410 --> 00:38:11.400
and the predicted travel
time in this [inaudible].

00:38:11.400 --> 00:38:17.350
And by minimizing this misfit function,
we like to update the velocity model.

00:38:17.350 --> 00:38:21.340
But, see, if you look at your data, you can 
see that a whole bunch of information

00:38:21.340 --> 00:38:25.200
that have not been taken into account.
And we know it’s very difficult to

00:38:25.200 --> 00:38:29.840
acquire those data, and we don’t want
to throw them away like garbage.

00:38:29.840 --> 00:38:32.560
So the full waveform inversion
is actually trying to take advantage

00:38:32.560 --> 00:38:36.860
of all of those energies.
So, in this case, instead of

00:38:36.860 --> 00:38:40.820
predict the travel times,
let’s predict the entire wavefield.

00:38:40.820 --> 00:38:44.890
So assuming that I have this very
smooth a priori velocity model,

00:38:44.890 --> 00:38:48.560
then we can solve the wave
equation through this velocity model

00:38:48.560 --> 00:38:51.640
and observe this
simulated wavefield.

00:38:51.640 --> 00:38:56.510
And then collect predicted data
at Earth surface, as shown here.

00:38:56.510 --> 00:39:00.750
So the goal of full waveform inversion
is actually to minimize the [inaudible]

00:39:00.750 --> 00:39:04.330
between your observed and
the predicted seismograms.

00:39:04.330 --> 00:39:09.410
So our goal is to actually minimize
observed data and predict the data

00:39:09.410 --> 00:39:14.140
by minimizing this least square function
and then find this optimal velocity

00:39:14.140 --> 00:39:17.920
model, which can minimize the
difference between these two.

00:39:17.920 --> 00:39:22.000
So next, I would like to test this
algorithm for 2D benchmark model.

00:39:22.010 --> 00:39:24.660
This is a 2D benchmark
Marmousi model.

00:39:24.660 --> 00:39:27.620
So this is the velocity
model we would like to recover.

00:39:27.620 --> 00:39:31.280
So, in which case,
you have this flat, sedimentary layer,

00:39:31.280 --> 00:39:35.820
and you have those faults.
And assuming that I have a relatively

00:39:35.820 --> 00:39:40.300
smooth a priori velocity model
as the starting model, as shown here,

00:39:40.300 --> 00:39:45.400
and if we run iteration for full waveform
inversion, you can see how this model

00:39:45.400 --> 00:39:50.740
evolves as we iteratively update it.
So you can see all these small details,

00:39:50.740 --> 00:39:55.230
they actually naturally emerge from this
smooth background velocity model.

00:39:55.230 --> 00:39:58.250
And, in this case, after
hundreds iterations, this is the

00:39:58.250 --> 00:40:02.620
velocity model we are able to build.
And if we compare this final model

00:40:02.620 --> 00:40:06.170
with the true model,
we’re able to see a lot of details.

00:40:06.170 --> 00:40:09.150
So on this – this point,
I just reach to 10 hertz.

00:40:09.150 --> 00:40:14.740
And if you want to see more details,
you can push to higher frequencies.

00:40:15.660 --> 00:40:20.830
So this is just [inaudible] data test.
We don’t have any noise.

00:40:20.830 --> 00:40:24.110
And the – really interesting
to draw the industry to pursue

00:40:24.110 --> 00:40:26.960
this technology is to
[inaudible] example.

00:40:26.960 --> 00:40:32.440
For instance, this is an example shown
by BP a couple of decades ago

00:40:32.440 --> 00:40:36.230
for the velocity model for
Valhall oil field in Norway.

00:40:36.230 --> 00:40:40.830
So those are the map view –
it’s 150 meter and 1,050 meter –

00:40:40.830 --> 00:40:43.600
by using classical
travel time tomography.

00:40:43.600 --> 00:40:46.840
So you can see the resolution
is actually not very good.

00:40:46.840 --> 00:40:50.120
However, once they tried full
waveform inversion, those are the

00:40:50.120 --> 00:40:55.610
velocity models that they were able to
build at shallow and greater depths.

00:40:55.610 --> 00:40:59.730
Well, in this case, we do observe
this channel field at shallow depths,

00:40:59.730 --> 00:41:05.130
and these – even this gas cloud
and even this fracture here,

00:41:05.130 --> 00:41:08.360
you can image by
using this technology.

00:41:09.160 --> 00:41:11.280
So we would like to use this
technology to build

00:41:11.280 --> 00:41:14.040
a crustal velocity
model for central Oklahoma.

00:41:14.040 --> 00:41:17.940
So at current stage, we are still at
the first iteration, and we compute

00:41:17.940 --> 00:41:21.350
those so-called sensitivity
kernel to see their behavior.

00:41:21.350 --> 00:41:26.610
For instance, this is how the sensitivity
kernel look like at 1-kilometer depth.

00:41:26.610 --> 00:41:31.550
And the color here – blue color
suggest that when we perform iteration,

00:41:31.550 --> 00:41:34.910
we should reduce the velocity.
And the red color suggest

00:41:34.910 --> 00:41:37.790
we should increase the
velocity in order to minimize

00:41:37.790 --> 00:41:41.180
a difference between
data and the prediction.

00:41:41.180 --> 00:41:49.140
Let’s see how it look like
at 3 kilometer, 5, 10, 15, and 20.

00:41:49.150 --> 00:41:53.880
So we can see, at first iteration,
we probably have very good constraint

00:41:53.880 --> 00:41:58.510
for the depths shallower than
10 to 15 kilometer, which is also

00:41:58.510 --> 00:42:01.720
the occurrence depths
for the seismicity in this area.

00:42:01.720 --> 00:42:05.000
So we would like to construct
a 3D velocity model with

00:42:05.000 --> 00:42:09.550
very good illumination at depths
shallower than 10 to 15 kilometer.

00:42:10.340 --> 00:42:14.660
So we did this before,
but this is for much larger scale.

00:42:14.670 --> 00:42:19.230
In this case, we did this for
north Texas and Oklahoma

00:42:19.230 --> 00:42:23.180
by using the same technology,
but we used ambient noise data.

00:42:23.180 --> 00:42:27.140
So, in this case, instead of
inverse dispersion relationship,

00:42:27.140 --> 00:42:30.330
we actually trying to invert
those empirical Green’s function

00:42:30.330 --> 00:42:35.700
by using full waveform inversion.
So those are the constructed

00:42:35.700 --> 00:42:42.770
tomographic model we have
at 3-, 6-, 10-, and 20-kilometer depths.

00:42:42.770 --> 00:42:50.030
And we also superimposed those
tectonic trajectory on top of these maps.

00:42:50.030 --> 00:42:53.140
And you can see the nice correlation
between our tomographic model

00:42:53.140 --> 00:42:57.440
with this geological structure.
For instance, if you compare it with

00:42:57.440 --> 00:43:03.240
tectonic map for Oklahoma and Texas,
for instance, we do observe this

00:43:03.240 --> 00:43:10.400
very sharp contrast across
this southern Oklahoma [inaudible]

00:43:10.400 --> 00:43:13.320
and Ouachita orogenic
front here.

00:43:13.320 --> 00:43:16.540
And if you look at this one,
you also see that our tomographic

00:43:16.540 --> 00:43:22.000
imaging can capture this line
of uplift in central Texas.

00:43:22.000 --> 00:43:24.800
And more interesting, if we
compare our seismic measurements

00:43:24.800 --> 00:43:29.050
with gravity data – for instance,
this is the gravity model we have.

00:43:29.050 --> 00:43:32.900
We do observe this correlation
between this fast velocity anomaly

00:43:32.900 --> 00:43:39.310
here with those gravity high
across Oklahoma and Texas.

00:43:39.310 --> 00:43:43.770
So this is very promising to show
that this technology really can help us

00:43:43.770 --> 00:43:49.160
to constrain or construct
3D crustal velocity models.

00:43:49.160 --> 00:43:54.500
So in conclusion, for elastic RTM,
we used the Helmholtz decomposition

00:43:54.500 --> 00:43:58.950
to separate P and S wavefields.
And also decomposed the P and S

00:43:58.950 --> 00:44:05.280
wavefield, the vector field, and have
correct amplitude as well as the phase.

00:44:05.280 --> 00:44:08.990
And then we introduced this vector
imaging condition to construct P-P

00:44:08.990 --> 00:44:12.990
and P-S images, and they
have high image quality

00:44:12.990 --> 00:44:17.540
and enable us to avoid
this polarity reversal issue.

00:44:18.040 --> 00:44:21.970
And also, we can use these separated
P and S waves to better locate

00:44:21.970 --> 00:44:28.280
earthquakes and imaging earthquake
rupture in 3D complex velocity models.

00:44:28.280 --> 00:44:32.670
And at current stage, we are trying to
apply this technology to real data.

00:44:32.670 --> 00:44:35.640
We are trying to apply
this for two sites – first,

00:44:35.640 --> 00:44:40.790
for southern California and
then for central Oklahoma.

00:44:40.790 --> 00:44:44.160
And we are going to use full
waveform inversion to map 3D

00:44:44.160 --> 00:44:47.200
elastic velocity distribution
in central Oklahoma,

00:44:47.200 --> 00:44:52.000
which might help us to better study
the seismicity in this region.

00:44:52.000 --> 00:44:55.120
And I would like to acknowledge
Texas Advanced Computing Center

00:44:55.120 --> 00:45:00.160
to provide computational resources
to perform those 3D calculations.

00:45:00.160 --> 00:45:02.040
And I would like to –
thanks for attention.

00:45:02.040 --> 00:45:05.440
I’d like to take questions
and comments. Thanks.

00:45:05.440 --> 00:45:11.360
[Applause]

00:45:12.540 --> 00:45:15.660
[Silence]

00:45:15.660 --> 00:45:17.040
Yeah.

00:45:17.960 --> 00:45:20.780
[Silence]

00:45:21.680 --> 00:45:25.890
- So first, just let me say, wow.
That was a lot. [laughter]

00:45:25.890 --> 00:45:29.000
A lot of different applications.
- This is just half, actually.

00:45:29.000 --> 00:45:31.300
- And great, great talk.
- Thanks.

00:45:31.300 --> 00:45:32.900
- I’m curious.
You know, you showed the

00:45:32.900 --> 00:45:36.060
example where you went to –
where you have 80% of the data.

00:45:36.070 --> 00:45:38.650
The reality probably is the
other way around, where we’re

00:45:38.650 --> 00:45:41.220
going to have about 20%.
- Right. Right.

00:45:41.220 --> 00:45:46.820
- And then also, the reality is, we have
a lot of uncertainty in our 3D models.

00:45:46.820 --> 00:45:55.320
So, you know, having a known synthetic
model that you’ve – you work with is

00:45:55.320 --> 00:46:00.010
one thing, but assuming we have
uncertainty in our velocity model …

00:46:00.010 --> 00:46:01.550
- Right.
- … and then only 20% of the data,

00:46:01.550 --> 00:46:04.690
how well do you think the
earthquake imaging and rupture imaging

00:46:04.690 --> 00:46:09.010
will do with these methods?
- Yeah. I didn’t show in this talk –

00:46:09.010 --> 00:46:13.430
actually, you know, I submitted a paper.
We did provide 20% – if you

00:46:13.430 --> 00:46:15.940
eliminate all – 80% data,
how does it look like?

00:46:15.940 --> 00:46:20.360
And the result is actually
very promising in that regard.

00:46:20.360 --> 00:46:25.040
And in terms of uncertainty for velocity,
at current stage, we don’t know.

00:46:25.040 --> 00:46:29.920
So unless we try this test for southern
California and central Oklahoma,

00:46:29.930 --> 00:46:34.480
then we will see how good we are.
At current stage, we – because

00:46:34.480 --> 00:46:37.220
right now, my grad students,
they are in Houston for internship.

00:46:37.220 --> 00:46:41.670
So once they come back, we can test
this for field data example. Yep.

00:46:41.670 --> 00:46:46.080
- Another thing you could do is tweak
your synthetic model to perturb it.

00:46:46.080 --> 00:46:48.260
- Sure. Yeah, yep, yep, yep.
- And then use that as the

00:46:48.260 --> 00:46:51.080
model that you try and …
- Yeah. We will try – say, for instance,

00:46:51.080 --> 00:46:59.030
add an error to your velocity models to
see how robust of our algorithm. Yeah.

00:46:59.030 --> 00:47:02.700
That is a good suggestion we will try.
Yep, thanks.

00:47:05.230 --> 00:47:07.800
- I echo Jamie. Wow,
that was really interesting.

00:47:07.800 --> 00:47:13.200
I had a question about the shallow
velocity model for central Oklahoma.

00:47:13.200 --> 00:47:18.160
There were some cool circles
in the model – if you go back

00:47:18.160 --> 00:47:21.620
a couple more slides.
It was the more detailed one.

00:47:21.620 --> 00:47:28.380
So the cool and hot colors are – do you
interpret those as wells – injection wells?

00:47:28.390 --> 00:47:30.920
Or do you – do you have an
interpretation of those?

00:47:30.920 --> 00:47:34.150
- Oh, no. Those are actually
the source and receiver locations.

00:47:34.150 --> 00:47:36.960
- Oh. Okay, okay.
- Yeah. Used to constrain

00:47:36.960 --> 00:47:39.200
those velocity models.
- Yeah, yeah. Okay. All right.

00:47:39.200 --> 00:47:43.880
- Yeah. So this is just for first iteration.
So once we – more interations,

00:47:43.880 --> 00:47:48.480
you will see better solution.
This one just tell us, how do you

00:47:48.480 --> 00:47:52.140
update your velocity
model for the first iteration.

00:47:52.140 --> 00:47:55.840
Suggest, for instance, as I mentioned,
blue color suggests that we

00:47:55.840 --> 00:47:58.740
need to reduce velocity.
And red color suggests we need to

00:47:58.740 --> 00:48:03.620
up the velocity in order to correct
those travel time differences. Yep.

00:48:03.620 --> 00:48:06.180
- All right. Thank you.
- Yep. Thanks.

00:48:11.480 --> 00:48:13.880
- That was a great talk.
Really interesting.

00:48:13.880 --> 00:48:15.880
- Thanks.
- And I was really struck by how

00:48:15.880 --> 00:48:18.490
well you could image that gas pocket.
So one thing I’m wondering,

00:48:18.490 --> 00:48:22.260
it’s been hard to image CO2
sequestration pilot projects,

00:48:22.260 --> 00:48:25.030
for example, when people
are injecting CO2 underground.

00:48:25.030 --> 00:48:28.050
Have you thought about using
full waveform migration and

00:48:28.050 --> 00:48:32.940
this technique to image the diffusion of
any gas – could be CO2 – anything else?

00:48:32.940 --> 00:48:37.720
- Sure, sure.
One faculty in our department,

00:48:37.720 --> 00:48:43.250
he is exploration seismologist.
He actually tried to use this technology

00:48:43.250 --> 00:48:46.080
to monitor – once you
inject the CO2,

00:48:46.080 --> 00:48:50.040
how the velocity perturbation
evolved as a function of time.

00:48:50.040 --> 00:48:52.850
So do this for 4D case
instead of 3D imaging.

00:48:52.850 --> 00:48:57.010
So they do – able to see
the evolution, migration,

00:48:57.010 --> 00:49:01.720
of those low-velocity heterogeneities,
yeah. People did try that.

00:49:01.720 --> 00:49:07.000
But in that case, is for, I think, for
[inaudible] tomography, so it’s 2D case.

00:49:07.000 --> 00:49:08.980
- Great.
- Yep.

00:49:08.980 --> 00:49:12.800
- And one other idea.
One place where we have a really good

00:49:12.800 --> 00:49:16.600
traditional 3D tomographic image is
at the SAFOD – San Andreas Fault

00:49:16.600 --> 00:49:19.720
Observatory at Depth site in Parkfield –
P and S wave tomography …

00:49:19.720 --> 00:49:22.510
- Yep, yep.
- … from 60 stations at the surface

00:49:22.510 --> 00:49:25.470
and very good images of
the rupture process for

00:49:25.470 --> 00:49:28.020
small earthquakes from
borehole seismometers.

00:49:28.020 --> 00:49:30.299
And so that would be something
really interesting to try.

00:49:30.300 --> 00:49:32.140
- Yeah, yeah, yeah. Yeah, yeah.
- So I’ll talk to you about it

00:49:32.140 --> 00:49:33.980
this afternoon.
- Yeah. Great, great, great. Thanks, thanks.

00:49:33.980 --> 00:49:37.200
Yeah, we would like to try that. Yeah.
- That would be great. Thanks.

00:49:40.080 --> 00:49:43.040
- Additional questions
before we break?

00:49:43.760 --> 00:49:45.360
Well, thank you.
That was a great talk.

00:49:45.360 --> 00:49:47.000
- Great. Yeah, thanks.
- And it’s a pleasure to have you.

00:49:47.000 --> 00:49:50.880
If you’d like to join us for lunch,
just come find us now-ish.

00:49:50.880 --> 00:49:53.220
And you’re all free to go.
Thank you.

00:49:53.660 --> 00:49:56.800
[Applause]

00:49:56.800 --> 00:50:03.700
[Silence]