Source code for grond.analysers.target_balancing.analyser

import copy
import time
import logging
import numpy as num
from pyrocko.guts import Int, Float, Bool
from pyrocko import gf

from grond.meta import Forbidden, GrondError

from ..base import Analyser, AnalyserConfig, AnalyserResult

logger = logging.getLogger('grond.analysers.target_balancer')


guts_prefix = 'grond'


[docs]class TargetBalancingAnalyser(Analyser): """ Estimating target weights that balance the signal amplitudes. Signal amplitudes depend on the source-receiver distance, on the phase type and the taper used. Large signals have in general a higher contribution to the misfit than smaller signals, without carrying more information. With this function, weights are estimated that shall balance the phase contributions. The weight estimation is based on synthetic waveforms that stem from a given number of random forward models. The inverse of the mean synthetic signal amplitudes gives the balancing weight. This is described as adaptive station weighting in Heimann (2011). """ def __init__(self, niter, use_reference_magnitude, cutoff): Analyser.__init__(self) self.niter = niter self.use_reference_magnitude = use_reference_magnitude self.cutoff = cutoff def log_progress(self, problem, iiter, niter): t = time.time() if self._tlog_last < t - 10. \ or iiter == 0 \ or iiter == niter - 1: logger.info( 'Target balancing for "%s" at %i/%i.' % ( problem.name, iiter, niter)) self._tlog_last = t def analyse(self, problem, ds): if self.niter == 0: return wtargets = [] if not problem.has_waveforms: return for target in problem.waveform_targets: wtarget = copy.copy(target) wtarget.flip_norm = True wtarget.weight = 1.0 wtargets.append(wtarget) wproblem = problem.copy() wproblem.targets = wtargets xbounds = wproblem.get_parameter_bounds() misfits = num.zeros((self.niter, wproblem.ntargets, 2)) rstate = self.get_rstate(problem) isbad_mask = None self._tlog_last = 0 for iiter in range(self.niter): self.log_progress(problem, iiter, self.niter) while True: if self.use_reference_magnitude: try: fixed_magnitude = wproblem.base_source.get_magnitude() except gf.DerivedMagnitudeError: raise GrondError( 'Cannot use use_reference_magnitude for this type ' 'of source model.') else: fixed_magnitude = None x = wproblem.random_uniform( xbounds, rstate, fixed_magnitude=fixed_magnitude) try: x = wproblem.preconstrain(x) break except Forbidden: pass if isbad_mask is not None and num.any(isbad_mask): isok_mask = num.logical_not(isbad_mask) else: isok_mask = None misfits[iiter, :, :] = wproblem.misfits( x, mask=isok_mask, nthreads=8) isbad_mask = num.isnan(misfits[iiter, :, 1]) mean_ms = num.mean(misfits[:, :, 0], axis=0) mean_ps = num.mean(misfits[:, :, 1], axis=0) weights = 1. / mean_ps families, nfamilies = wproblem.get_family_mask() for ifamily in range(nfamilies): weights[families == ifamily] /= ( num.nansum(weights[families == ifamily]) / num.nansum(num.isfinite(weights[families == ifamily]))) if self.cutoff is not None: weights[mean_ms / mean_ps > self.cutoff] = 0.0 for weight, target in zip(weights, problem.waveform_targets): target.analyser_results['target_balancing'] = \ TargetBalancingAnalyserResult(weight=float(weight)) for itarget, target in enumerate(problem.waveform_targets): logger.info(( 'Balancing analysis for target "%s":\n' ' m/p: %g\n' ' weight: %g\n' ) % ( target.string_id(), mean_ms[itarget] / mean_ps[itarget], weights[itarget]))
class TargetBalancingAnalyserResult(AnalyserResult): weight = Float.T()
[docs]class TargetBalancingAnalyserConfig(AnalyserConfig): """Configuration parameters of the target balancing.""" niterations = Int.T( default=1000, help='Number of random forward models for mean phase amplitude ' 'estimation') use_reference_magnitude = Bool.T( default=False, help='Fix magnitude of random sources to the magnitude of the ' 'reference event.') cutoff = Float.T( optional=True, help='Remove targets where ratio m/p > cutoff, where m is the misfit ' 'between synthetics and observations and p is the misfit between ' 'synthetics and zero-traces. Magnitude should be fixed to use ' 'this.') def get_analyser(self): return TargetBalancingAnalyser( niter=self.niterations, use_reference_magnitude=self.use_reference_magnitude, cutoff=self.cutoff)
__all__ = ''' TargetBalancingAnalyser TargetBalancingAnalyserConfig '''.split()