Earthquake location by backprojection

What would happen if time could run backwards as well as forwards ? 

Well seismologists can run time backwards in a computer and work out where an earthquake has happened. 


We know what time a seismic wave reaches a monitoring station, and we know the speed it travels.  If we run time backwards we can then work out the locus of possible origins for the source of the earthquake at each time in the past, for each station this will be a circle of possible origin locations that increases in size as we go back in time.   With three or more stations the origin location (and time) of the event can be estimated from the point when all of these possible locations converge on one point.

A simple python program and some data from a recent earthquake in Cornwall illustrates this principle
import matplotlib.pyplot as plt
import math
import numpy as np
from matplotlib.patches import Circle
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation

#Writer = animation.writers['ffmpeg']
#writer = Writer(fps=20, metadata=dict(artist='Paul Denton'),bitrate=1800)

#image_file = background image of cornwall
# background plot 20px/km
scale=20
img = plt.imread('./corn_map.png')
stations=['PZ','CCA1','RL','FL']
x=[315,764,946,958]
y=[487,330,164,417]
#pick times (sec after 16:52
pick=[8.430,7.600,9.418,7.880]
# Create a figure. Equal aspect so circles look circular
fig,ax = plt.subplots(1)
ax.set_aspect('equal')
# pwave velocity in crust =5.76km/sec
speed=5.76
#depth=10km
depth=10.0
#set latest possible time for event (directly below closest station)
latest_time=min(pick)-depth/speed -0.01
# Create a figure. Equal aspect so circles look circular
fig,ax = plt.subplots(1)
ax.set_aspect('equal')
ax.imshow(img)


def run_animation():
    anim_running = True
    # section to allow stopping of animation ona key press
    def onClick(event):
        nonlocal anim_running
        if anim_running:
            anim.event_source.stop()
            anim_running = False
        else:
            anim.event_source.start()
            anim_running = True

    def update(frame_number):

    #assume 1x fframenumber=10 msec
        ax.clear()
        ax.imshow(img)
        #  set a test origin time and go backwards in time
        time=round(latest_time-(frame_number/100),2)
        label='time=16:52:'+str(time)
        # plot time label on animation
        ax.text(100,750, label)
        for i in range(4):
            traveltime=pick[i]-time  # travel time to each station for this test origin
            slope=traveltime*speed   #  slopw distance from station for this test origin
            epidist=math.sqrt(slope*slope-depth*depth)   #   surface distance for this slope distance
            size=epidist*scale # convert surface distance to map scale
            circ = Circle((x[i],y[i]),size, color='k',fill=False)
            # plot a range sircle of all possible locations that match
            # this test origin time for each station
            ax.add_patch(circ)

    fig.canvas.mpl_connect('button_press_event', onClick)  # pause animation when a mouse is clicked

    anim = animation.FuncAnimation(fig, update,frames=200)

run_animation()

# Show the image
#animation.save('cornwall_backprojection.mp4')
plt.show()

more information and resources on educational and citizen seismology at https://www.dentonseismo.co.uk 

Comments

Popular posts from this blog

lego seismometer

Getting data from a BBC:microbit

Barnsley Tyke quakes