Plot objective function

The output file objective_function.txt contains important information, which directly reflects the performance of inversion. Here, we demonstrate how to visualize the objective function reduction curve and how to read the underlying information.

Here, we demonstrate how to plot earthquakes and stations. First, we load necessary Python modules


import sys
sys.path.append('../utils')
import functions_for_data as ffd
 
import os
try: 
    os.mkdir('img')
except:
    pass

1. Read file

We first the file OUTPUT_FILES/objective_function.txt.


path = "OUTPUT_FILES"
full_curve, location_curve, model_curve = ffd.read_objective_function_file(path)
 
print("full_curve: ", full_curve.shape, ", the total objective function value during the inversion, including relocation and model update")
print("location_curve: ", location_curve.shape, ", the objective function value during the relocation step")
print("model_curve: ", model_curve.shape, ", the objective function value during the model update step")
 
print("The first index is iteraion number, the second index is the objective function value vector")

full_curve contains the information at every step;
location_curve only contains the information with the type of relocation;
model_curve only contains the information with the type of model update;

2. Information in the objective function file

(a). The objective function value

We can get the objective function value used for inversion at every iterations by


# (Option 1) objective function value
full_obj = full_curve[:,0]  
location_obj = location_curve[:,0]
model_obj = model_curve[:,0]

The objective function value directly reflects the overall consistency between observed and syntehtic traveltime data used for model update and/or relocation.
Depening on the used data type, the objective function value may be contributed by one or several datasets, including the absolute traveltimes, common-source differential arrival times, common-receiver differential arrival times, and teleseismic relative traveltimes.
Typically, in real-data inversion, the objective function value keeps decreasing in the beginning sevel or tens iterations, and then steadily converges to a specific value due to intrinsic generalized noise.
In synthetic tests, the objective function value may decrease to the noise level if the input anomalies are largely resolved.

(b) The objective function value of single data type

We can extract the objective function value of absolute traveltimes:


# (Option 2) objective function value for only traveltime
full_obj_tt = full_curve[:,1]
location_obj_tt = location_curve[:,1]
model_obj_tt = model_curve[:,1]

objective function value of common-source differential arrival times


# (Option 3) objective function value for only common source differential arrival time
full_obj_cs = full_curve[:,2]
location_obj_cs = location_curve[:,2]
model_obj_cs = model_curve[:,2]

objective function value of common-receiver differential arrival times


# (Option 4) objective function value for only common receiver differential arrival time
full_obj_cr = full_curve[:,3]
location_obj_cr = location_curve[:,3]
model_obj_cr = model_curve[:,3]

objective function value of teleseismic differential arrival time


# (Option 5) objective function value for teleseismic differential arrival time
full_obj_tele = full_curve[:,4]
location_obj_tele = location_curve[:,4]
model_obj_tele = model_curve[:,4]

Specific dataset may be included or excluded in model update and/or relocation based on the settings in the YAML file.


###################################################
#          model update parameters setting        #
###################################################
model_update:  
  ...
# -------------- using absolute traveltime data --------------
  abs_time:
    use_abs_time: true # 'true' for using absolute traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
 
  # -------------- using common source differential traveltime data --------------
  cs_dif_time:
    use_cs_time: false # 'true' for using common source differential traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
 
  # -------------- using common receiver differential traveltime data --------------
  cr_dif_time:
    use_cr_time: false # 'true' for using common receiver differential traveltime data to update model parameters; 'false' for not using (no need to set parameters in this section)
  ...
 
#################################################
#          relocation parameters setting        #
#################################################
relocation: # update earthquake hypocenter and origin time (when run_mode : 2 and 3)
  ...
  # -------------- using absolute traveltime data --------------
  abs_time:
    use_abs_time : false # 'yes' for using absolute traveltime data to update ortime and location; 'no' for not using (no need to set parameters in this section)
 
  # -------------- using common source differential traveltime data --------------
  cs_dif_time:
    use_cs_time : false # 'yes' for using common source differential traveltime data to update ortime and location; 'no' for not using (no need to set parameters in this section)
 
  # -------------- using common receiver differential traveltime data --------------
  cr_dif_time:
    use_cr_time : false # 'yes' for using common receiver differential traveltime data to update ortime and location; 'no' for not using (no need to set parameters in this section)

(c) Mean value of data residuals

Data residual means the difference between synthetic and observed data.

We can get the mean value of all data residuals at every iterations by


# (Option 6) mean value of all data residual
full_mean = full_curve[:,5]
location_mean = location_curve[:,5]
model_mean = model_curve[:,5]

Similarly, mean value of data residuals for single data type can be also obatined

for absolute traveltimes


# (Option 8) mean value of residuals of traveltime
full_mean_tt = full_curve[:,7]
location_mean_tt = location_curve[:,7]
model_mean_tt = model_curve[:,7]

The mean value of data residuals for absolute traveltime is an important indicator. A positive value means that the synthetic traveltimes are systematically larger than the observed traveltimes, indicating that the model velocity is too low under the earthquake locations, or the facal depths are too deep under the velocity model, and vice versa.

Typically, it is important (or necessary) to guarantee the mean value of absolute traveltime residuals is near 0 in the beginning of the inversion, e.g., emperically ranging from -0.15 to 0.15. Otherwise, the model will be systematically deviated upward or downward in the following iterations, increasing the risk of being trapped in local minima.

for common-source differential arrival times


# (Option 10) mean value of residuals of common source differential arrival time
full_mean_cs = full_curve[:,9]
location_mean_cs = location_curve[:,9]
model_mean_cs = model_curve[:,9]

for common-receiver differential arrival times


# (Option 12) mean value of residuals of common receiver differential arrival time
full_mean_cr = full_curve[:,11]
location_mean_cr = location_curve[:,11]
model_mean_cr = model_curve[:,11]

for teleseismic differential arrival times


# (Option 14) mean value of residuals of teleseismic differential arrival time
full_mean_tele = full_curve[:,13]
location_mean_tele = location_curve[:,13]
model_mean_tele = model_curve[:,13]

(d) Standard deviation of data residuals

A simple but rough understanding of the standard deviation of data residuals is the overall difference in second between the synthetic and observed traveltimes.

Typically, the standard deviation is similar to the squared root of the objection function if a consistent weight is used for all data.

We can get the standard deviation of all data residuals at every iterations by


# (Option 7) standard deviation of all data residual
full_std = full_curve[:,6]
location_std = location_curve[:,6]
model_std = model_curve[:,6]

and get the standard deviation of single data type

for absolute traveltimes


# (Option 9) standard deviation of residuals of traveltime
full_std_tt = full_curve[:,8]
location_std_tt = location_curve[:,8]
model_std_tt = model_curve[:,8]

for common-source differential arrival times


# (Option 11) standard deviation of residuals of common source differential arrival time
full_std_cs = full_curve[:,10]
location_std_cs = location_curve[:,10]
model_std_cs = model_curve[:,10]

for common-receiver differential arrival times


# (Option 13) standard deviation of residuals of common receiver differential arrival time
full_std_cr = full_curve[:,12]
location_std_cr = location_curve[:,12]
model_std_cr = model_curve[:,12]

for teleseismic differential arrival times


# (Option 15) standard deviation of residuals of teleseismic differential arrival time
full_std_tele = full_curve[:,14]
location_std_tele = location_curve[:,14]
model_std_tele = model_curve[:,14]

1. Plot objective function reduction curve

For example, we plot the objective function reduction curve related to model update.


# %%
# plot objective functin reduction
 
import matplotlib.pyplot as plt
import numpy as np
 
plt.figure(figsize=(10, 6))
ax = plt.subplot(1, 1, 1)
 
ax.plot(model_obj/np.max(model_obj), label='objective function', linewidth=2)
ax.set_xlim([-0.2, len(model_obj)-0.8])
ax.set_ylim([0, 1.1])
ax.grid()
ax.set_xlabel('Iteration number',fontsize=14)
ax.set_ylabel('Normalized value',fontsize=14)
ax.tick_params(axis='x', labelsize=14)  
ax.tick_params(axis='y', labelsize=14) 
ax.legend(fontsize=14)
 
plt.savefig('img/5_objective_function_reduction.png', dpi=300, bbox_inches='tight', edgecolor='w', facecolor='w')