///|
/// Configuration for a weighted hypocenter search.
pub struct LocationSearchConfig {
latitude_min : Double
latitude_max : Double
longitude_min : Double
longitude_max : Double
depth_min_km : Double
depth_max_km : Double
grid_size : Int
depth_levels : Int
maximum_residual : Double
}
///|
/// Residual statistics for a location solution.
pub struct LocationResidualSummary {
rms : Double
median_absolute : Double
maximum_absolute : Double
inlier_count : Int
outlier_count : Int
score : Double
}
///|
/// Approximate uncertainty ellipse around a location.
pub struct LocationUncertainty {
horizontal_sigma_km : Double
latitude_sigma_degrees : Double
longitude_sigma_degrees : Double
depth_sigma_km : Double
azimuth_degrees : Double
confidence_level : Double
}
///|
/// A location result with residual and uncertainty diagnostics.
pub struct RefinedLocation {
estimate : LocationEstimate
residuals : LocationResidualSummary
uncertainty : LocationUncertainty
iterations : Int
converged : Bool
}
///|
/// Return a conservative search configuration around a source guess.
pub fn default_location_search(
latitude : Double,
longitude : Double,
) -> LocationSearchConfig {
{
latitude_min: latitude - 0.5,
latitude_max: latitude + 0.5,
longitude_min: longitude - 0.5,
longitude_max: longitude + 0.5,
depth_min_km: 0.0,
depth_max_km: 50.0,
grid_size: 9,
depth_levels: 6,
maximum_residual: 5.0,
}
}
///|
/// Compute weighted residual summary with a Huber-like outlier policy.
pub fn summarize_location_residuals(
residuals : Array[Double],
maximum_residual : Double,
) -> LocationResidualSummary {
let limit = maximum_residual.abs()
let inliers = residuals.filter(value => value.abs() <= limit)
let outliers = residuals.length() - inliers.length()
let absolute = residuals.map(value => value.abs())
let score = if residuals.is_empty() {
0.0
} else {
inliers.length().to_double() /
residuals.length().to_double() *
@math.exp(-rms_samples(residuals) / limit.max(1.0e-12))
}
{
rms: rms_samples(residuals),
median_absolute: median_samples(absolute),
maximum_absolute: maximum_sample(absolute),
inlier_count: inliers.length(),
outlier_count: outliers,
score,
}
}
///|
/// Estimate the origin time by minimizing a robust residual objective.
pub fn robust_origin_time(
observed : Array[Double],
predicted_without_origin : Array[Double],
iterations : Int,
) -> Double {
let count = observed.length().min(predicted_without_origin.length())
if count == 0 {
0.0
} else {
let mut origin = median_samples(observed) -
median_samples(predicted_without_origin)
let rounds = iterations.max(1)
for _ in 0.. value.abs())).max(
1.0e-6,
)
let mut correction = 0.0
let mut weight = 0.0
for residual in residuals {
let w = (scale / residual.abs().max(scale)).min(1.0)
correction += residual * w
weight += w
}
if weight > 0.0 {
origin += correction / weight
}
}
origin
}
}
///|
/// Search latitude, longitude, and depth with a coarse-to-fine grid.
pub fn grid_search_hypocenter(
model : VelocityModel,
stations : Array[StationPosition],
observed : Array[Double],
config : LocationSearchConfig,
) -> RefinedLocation {
let grid = config.grid_size.max(2)
let depths = standard_depth_samples(
config.depth_max_km.max(config.depth_min_km),
config.depth_levels.max(1),
).map(value => value + config.depth_min_km)
let mut best = {
latitude: config.latitude_min,
longitude: config.longitude_min,
depth_km: config.depth_min_km,
origin_time: minimum_sample(observed),
residual: 1.0e300,
}
for depth in depths {
for lat_index = 0; lat_index < grid; lat_index = lat_index + 1 {
let latitude = config.latitude_min +
(config.latitude_max - config.latitude_min) *
lat_index.to_double() /
(grid - 1).to_double()
for lon_index = 0; lon_index < grid; lon_index = lon_index + 1 {
let longitude = config.longitude_min +
(config.longitude_max - config.longitude_min) *
lon_index.to_double() /
(grid - 1).to_double()
let source = StationPosition::new("source", latitude, longitude, 0.0)
let predicted = Array::new()
for station_index = 0
station_index < stations.length() &&
station_index < observed.length()
station_index = station_index + 1 {
let distance = surface_distance_km(source, stations[station_index])
predicted.push(travel_time(model, depth, distance, false))
}
let origin = robust_origin_time(observed, predicted, 3)
let residual = location_residual(
model, depth, latitude, longitude, origin, stations, observed,
)
if residual < best.residual {
best = {
latitude,
longitude,
depth_km: depth,
origin_time: origin,
residual,
}
}
}
}
}
let source = StationPosition::new(
"source",
best.latitude,
best.longitude,
0.0,
)
let predicted = Array::new()
for station in stations {
let distance = surface_distance_km(source, station)
predicted.push(
best.origin_time + travel_time(model, best.depth_km, distance, false),
)
}
let residuals = Array::new()
for index = 0
index < observed.length() && index < predicted.length()
index = index + 1 {
residuals.push(observed[index] - predicted[index])
}
let summary = summarize_location_residuals(residuals, config.maximum_residual)
let horizontal = (config.latitude_max - config.latitude_min).abs() *
111.32 /
grid.to_double().max(1.0)
let uncertainty = {
horizontal_sigma_km: horizontal,
latitude_sigma_degrees: horizontal / 111.32,
longitude_sigma_degrees: horizontal / 111.32,
depth_sigma_km: (config.depth_max_km - config.depth_min_km).abs() /
config.depth_levels.max(1).to_double(),
azimuth_degrees: 0.0,
confidence_level: summary.score,
}
{
estimate: best,
residuals: summary,
uncertainty,
iterations: grid * grid * depths.length(),
converged: summary.score >= 0.5,
}
}
///|
/// Refine a location by searching a shrinking neighborhood around an estimate.
pub fn refine_location_neighborhood(
model : VelocityModel,
stations : Array[StationPosition],
observed : Array[Double],
initial : LocationEstimate,
latitude_step : Double,
longitude_step : Double,
depth_step : Double,
iterations : Int,
) -> RefinedLocation {
let mut best = initial
let mut lat_step = latitude_step.abs()
let mut lon_step = longitude_step.abs()
let mut dep_step = depth_step.abs()
let mut converged = false
let rounds = iterations.max(1)
for _ in 0.. String {
"lat=" +
location.estimate.latitude.to_string() +
" lon=" +
location.estimate.longitude.to_string() +
" depth=" +
location.estimate.depth_km.to_string() +
" origin=" +
location.estimate.origin_time.to_string() +
" residual=" +
location.estimate.residual.to_string() +
" score=" +
location.residuals.score.to_string() +
" sigma_km=" +
location.uncertainty.horizontal_sigma_km.to_string() +
" converged=" +
location.converged.to_string()
}