Add option for individual travel time tables per station

[yetinam]

Currently, PyOcto uses the same travel time table for each stations. However, technically it would be pretty straightforward to use individual tables. This is useful to properly model station elevations. The current approximation for station elevation correction is inexact. While this is okay for classical tectonic settings, the approximation becomes problematic when stations and events are at similar depth, e.g., for induced seismicity or OBS sensors.

[zhong-yy]

Hi, using individual travel times is an interesting feature. I have followed this issue since it's created. Is there any news?

I think there may be different ways to implement this.

(1) The simplest way may be shifting the z coordinate of the volume in the beginning of the bool VelocityModel1D::contains(const Volume &volume, const Pick *pick) function and the double VelocityModel1D::travel_time(const Volume &volume, const std::string &station, char phase)

pyocto/src/VelocityModel.cpp

Line 183 in 32379f3

bool VelocityModel1D::contains(const Volume &volume, const Pick *pick) {

For example (I am not sure),

double VelocityModel1D::travel_time(const Volume &volume, const std::string &station, char phase){
  auto fake_volume=Volumne(volume);
  // use abs() to swap the source and receiver when the source is above the station
  fake_volume.z=std::abs(volume.z-station.z); 
  ...
}

But this implicitly assumes different 1D layered models for different stations because of different station elevations.

(2) Saving individual travel time tables for each station may be a bit more tricky and requires more changes to the code, but is more reasonable than the above option. To save individual time tables, as far as I understand, the p_times, s_times, local_minima_p, local_minima_s, local_maxima_p, local_maxima_s in the VelocityModel1D class in VelocityModel.h/cpp

pyocto/src/VelocityModel.h

Lines 88 to 89 in 32379f3

	double *p_times; // 2D array of travel times
	double *s_times; // 2D array of travel times

need to be changed to a map that takes station name as the key.

In addition, a reference depth should be specifed as the top surface of the model (currently, 0 is assumed as the top).

In the VelocityModel1D.create_model(...) method in associator.py, the source is assumed at 0 km depth,

pyocto/src/pyocto/associator.py

Line 164 in 32379f3

p_times[0, 0] = 0

It should be changed according to the station elevation and the reference depth of the the model top surface, for example (not sure if it makes sense)

p_times[int((station.z-reference_depth)/delta), 0] = 0

Meanwhile, the elevation correction in

pyocto/src/VelocityModel.cpp

Line 178 in 32379f3

return -station.z / surface_p_velocity;

should become

int sta_z_id=int((station.z-reference_depth)/delta);
double correction = (sta_z_id*delta-station.z)/surface_p_velocity;

(3) Maybe another way to do this is to use travel time tables in NonLinLoc. This may be the best solution. In this way, the same travel time tables could be used for association (pyocto), location (NonLinLoc) and relocation (GrowClust)

By the way, I have a question about the code of checking whether a volume contains a pick. In the following code, why do you compare an integer index z and the actual depth range volume.z_min() and volume.z_max()?

pyocto/src/VelocityModel.cpp

Lines 241 to 248 in 32379f3

	for (auto z : local_minima) {
	if (volume.z_min() <= z and z <= volume.z_max())
	min_candidates.push_back(times[index(ix_min, z)]);
	}
	for (auto z : local_maxima) {
	if (volume.z_min() <= z and z <= volume.z_max())
	max_candidates.push_back(times[index(ix_max, z)]);
	}

[yetinam]

Hey @zhong-yy ,

first of all, thanks a lot for spotting the bug with the 1D model that has been hidden in PyOcto for almost 3 years. That was indeed a serious issue. I've fixed it now in #31.

Regarding the individual models for each station, I unfortunately did not have time do work on this topic, but I'd be happy if you want to pick it up and submit a PR. The solution I'd recommend is (2), because this is the clean solution. In principle, there's nothing wrong about (3), but it would still only work for 1D tables and the overhead for reading NonLinLoc tables might be unreasonably large. 3D tables don't work, because PyOcto needs to be able to efficiently query what the minimum and maximum travel time from a volume to a station is. For all I know, there's no efficient solution for this for general 3D velocity models (but I'd be happy to be proven wrong).

For (2), what I would do is using a different "virtual" z=0 for each station, i.e., define for each station that it is located at z=0. Then, you can calculate a travel time table for the station by simply shifting all the layers in the velocity model down by z before running the eikonal solver. During association, you just need to transform the true z coordinate of the volume to the one with reference to the station depth. In any case, I think it's best if the station is located exactly on a grid point (by shifting the grid appropriately), because otherwise there'll always be unnecessary fractional errors.

A bit of care needs to be taken to include the possibility of events above stations. This matters, for example, when using borehole stations. The easiest solution is probably to add a constant number of grid nodes above the station. That shouldn't be difficult from a technical perspective but it might need some careful checking to make sure it's correctly added everywhere.

To keep downward compatibility, one option might be to implement a new class like StationSpecificVelocityModel1D and keep the existing classes around.

[zhong-yy]

Hello, @yetinam , thank you for the fix!

It is a good idea to implement a new class StationSpecificVelocityModel1D. I would be happy to submit a PR.

Then, you can calculate a travel time table for the station by simply shifting all the layers in the velocity model down by z before running the eikonal solver.

For shifting the layers, do you mean extending the first velocity layer up to the station position if the station is above the velocity model and removing the layers above the station if it is embedded in the velocity model?

Currently, if zlim is less than 0, volumes above 0 would be shifted to 0 in calculation:

pyocto/src/VelocityModel.cpp

Lines 227 to 228 in 48d4268

	if (iz_min < 0)
	iz_min = 0;

In the new class, we can remove this limit.

To add a constant number of grid notes, maybe we can add a parameter z_padding to the create_model() method, and save the number of nodes in the class. The number will be added to the z index when the travel time table is looked up.

It seems that there is another minor bug in the following code of creating velocity grids:

pyocto/src/pyocto/associator.py

Lines 175 to 180 in 48d4268

	p_speeds[:, int(depth1 / delta) : int(depth2 / delta)] = np.linspace(
	vp1, vp2, int(depth2 / delta) - int(depth1 / delta)
	)
	s_speeds[:, int(depth1 / delta) : int(depth2 / delta)] = np.linspace(
	vs1, vs2, int(depth2 / delta) - int(depth1 / delta)
	)

By default, numpy includes the end point in the resulting array. Here, vp2 and vs2 are assigned to one grid node before depth2, which causes a step at depth2. Maybe endpoint=False should have been used.

[yetinam]

For shifting the layers, do you mean extending the first velocity layer up to the station position if the station is above the velocity model and removing the layers above the station if it is embedded in the velocity model?

Exactly, that's what I would propose. I don't think you even need to remove the layers above if you allow for events above the station (see below).

To add a constant number of grid notes, maybe we can add a parameter z_padding to the create_model() method, and save the number of nodes in the class. The number will be added to the z index when the travel time table is looked up.

I was thinking of something like introducing z_padding_layers, such that the p_times[0, z_padding_layers] = 0.0. This means, there will always be delta * z_padding_layers kilometers above each station available in the search space.

Thanks also for pointing out the issue with the interpolation. Indeed, it should include endpoint=False. I'll fix that in a commit now. Fortunately, the difference shouldn't be huge, so I don't think it'll cause major issues.

[zhong-yy]

Sounds good! I will find some time to give it a try, but I would be inclined to make z_padding_layers have the same unit (i.e. km) as zdist, and calculate the number of padding layers the same way as calculating nx and nz:

pyocto/src/pyocto/associator.py

Lines 152 to 153 in 48d4268

	nx = int(xdist / delta)
	nz = int(zdist / delta)

Do you have any idea of how to test the new feature?

[zhong-yy]

Hello @yetinam , I have created a PR, could you review my code? I used the same example as the one for VelocityModel1D, but maybe this is not a good example for StationSpecificVelocityModel1D because most events are deep.

When I test it using real picks from a volcano, the initial locations of the associated events become more compact, but also become shallower. I am not sure whether it is a bug. I have checked the code carefully, but did not find a bug. In additiion, StationSpecificVelocityModel1D seems much slower than VelocityModel1D.

BTW, your code is really well written and easy to read :)

[zhong-yy]

Hello @yetinam , I have created a PR, could you review my code? I used the same example as the one for VelocityModel1D, but maybe this is not a good example for StationSpecificVelocityModel1D because most events are deep.

When I test it using real picks from a volcano, the initial locations of the associated events become more compact, but also become shallower. I am not sure whether it is a bug. I have checked the code carefully, but did not find a bug. In additiion, StationSpecificVelocityModel1D seems much slower than VelocityModel1D.

BTW, your code is really well written and easy to read :)

The shallower locations were caused by incorrect elevation correction. After I disabled elevation correction, the result became normal.