source: trunk/src/seismic_processing/glass/src/modules/OPCalc/opcalc.cpp @ 3212

Revision 3212, 16.7 KB checked in by paulf, 12 years ago (diff)

sync from hydra_proj circa Oct 2007

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
RevLine 
[2176]1#include <math.h>
2#include <Debug.h>
3#include <string.h>
4
5
6# define _OPCALC_SO_EXPORT  __declspec( dllexport)
7
8#include <opcalc.h>
9#include <opcalc_const.h>
10
[3212]11extern "C" {
12void __cdecl qsortd (double *base, unsigned num);
13}
14
15
16int OPCalc_ProcessOriginSimple(ORIGIN * pOrg, PICK ** pPick, int nPck);
17
18
[2176]19// Lib Global variables - File Scope
20static ITravelTime * pTT = 0x00;      // file-scope master traveltime pointer
21
22//---------------------------------------------------------------------------------------Zerg
23// Calculate the OPCalc_BellCurve function of x.
24// The height of the curve ranges from 0 to 1.
25// The curve is centered at 0, with a range from -1 to 1,
26// such that OPCalc_BellCurve(0) = 1 and OPCalc_BellCurve(-1) = OPCalc_BellCurve(1) = 0.
27// OPCalc_BellCurve is usually called as OPCalc_BellCurve(x/w),
28// where x is the deviation from the mean, and w is the half-width
29// of the distrubution desired.
30// OPCalc_BellCurve(-w) = OPCalc_BellCurve(w) = 0.0.
31// dBellCurve/dX(0) = dBellCurve/dX(1) = 0.0.
32// OPCalc_BellCurve is symetric about 0.0, and is roughly gaussian with a
33// half-width of 1. For x < -1 and x > 1, OPCalc_BellCurve is defined as 0.
34// In Carl terms, OPCalc_BellCurve is defined as Zerg.
35double OPCalc_BellCurve(double xx) 
36{
37        double x = xx;
38        if(x < 0.0)
39                x = -x;
40        if(x > 1.0)
41                return 0.0;
42        return 2.0*x*x*x - 3.0*x*x + 1.0;
43}
44
45
[3212]46int __cdecl ComparePickAffinityAsc(const void * pvoid1, const void * pvoid2)
47{
48  PICK * p1, *p2;
49
50  p1 = *(PICK **) pvoid1;
51  p2 = *(PICK **) pvoid2;
52
53  if(p1->dAff < p2->dAff)
54    return(-1);
55  if(p1->dAff > p2->dAff)
56    return(1);
57  else
58    return(0);
59}
60
61int __cdecl ComparePickResidDesc(const void * pvoid1, const void * pvoid2)
62{
63  PICK * p1, *p2;
64
65  p1 = *(PICK **) pvoid1;
66  p2 = *(PICK **) pvoid2;
67
68  if(fabs(p1->tRes) > fabs(p2->tRes))
69    return(-1);
70  if(fabs(p1->tRes) < fabs(p2->tRes))
71    return(1);
72  else
73    return(0);
74}
75 
76int __cdecl ComparePickiPickAsc(const void * pvoid1, const void * pvoid2)
77{
78  PICK * p1, *p2;
79
80  p1 = *(PICK **) pvoid1;
81  p2 = *(PICK **) pvoid2;
82
83  if(p1->iPick < p2->iPick)
84    return(-1);
85  if(p1->iPick > p2->iPick)
86    return(1);
87  else
88    return(0);
89}
90 
91
92int OPCalc_ProcessOrigin(ORIGIN * pOrg, PICK ** pPick, int nPck) 
93{
94
95  double dFullAff;
96  int    rc;
97  int    i;
98  PICK * pck;
99
100  /* process the origin with the full complement of picks
101   *******************************************************/
102  if(rc=OPCalc_ProcessOriginSimple(pOrg, pPick, nPck))
103    return(rc);
104
105  dFullAff = pOrg->dAff;
106
107  // log the results
108  CDebug::Log(DEBUG_MINOR_INFO,"OPC_PO: Full Origin: %.2f/%.2f/%.2f - %.2f  Aff: %.2f  AfGp: %.2f AfNeq: %.2f  RMS: %.2f Gap: %.2f\n",
109              pOrg->dLat, pOrg->dLon, pOrg->dZ, pOrg->dT, pOrg->dAff, pOrg->dAffGap, pOrg->dAffNumArr, pOrg->dRms, pOrg->dGap,
110              pOrg->dRms, pOrg->nEq);
111
112  for(i=0; i<nPck; i++)
113  {
114                pck = pPick[i];
115    CDebug::Log(DEBUG_MINOR_INFO,
116                "%3d) %-4s %-6s %5.1f %3.0f %3.0f %5.2f %6.1f %4.1f %3.1f/%3.1f/%3.1f  %s\n",
117                                  i, pck->sSite, pck->sPhase,
118                                  pck->dDelta, pck->dAzm, pck->dToa * RAD2DEG, pck->tRes, pck->dT - pOrg->dT, 
119                pck->dAff, pck->dAffDis, pck->dAffRes, pck->dAffPPD, pck->sTag
120               ); 
121  }
122
123  // sort the picks by resid
124  qsort(pPick, nPck, sizeof(PICK *), ComparePickResidDesc);
125
126  /* when trimming picks, don't trim anything that is either:
127     1) outside of the worst 20% in terms of residual
128     2) inside 1 std. deviation (assuming normal distribution with stddev = 3.33 sec)
129     We Make three BIG assumptions here:
130     a) That the ORIGIN is an accurate approximation of the Event origin
131     b) That the pick residual distribution is NORMAL.
132     c) That the stddev for such a distribution is 3.33 seconds.
133   **************************************************************/
134#define OPCALC_MAX_RESID_TRIM_RATIO .20
135#define OPCALC_PICK_RESID_DIST_STD_DEV_SEC 3.33
136
137  // trim out the worst picks (by res)
138  for(i=0; i<(nPck*OPCALC_MAX_RESID_TRIM_RATIO); i++)
139  {
140                pck = pPick[i];
141    CDebug::Log(DEBUG_MINOR_INFO, "Pick filter: %d res: %.2f aff: %.2f  (%.2f ~ %.2f)\n",
142                i, pck->tRes, pck->dAff,
143                fabs(pck->tRes), OPCALC_PICK_RESID_DIST_STD_DEV_SEC);
144    if(fabs(pck->tRes) < OPCALC_PICK_RESID_DIST_STD_DEV_SEC)
145      break;
146  }
147
148  /* process the origin with the full complement of picks
149   *******************************************************/
150  if(i<nPck)
151  {
152    OPCalc_ProcessOriginSimple(pOrg, &pPick[i], nPck-i);
153
154    // log the results
155    CDebug::Log(DEBUG_MINOR_INFO,"OPC_PO: Trimmed Origin: %.2f/%.2f/%.2f - %.2f  Aff: %.2f  AfGp: %.2f AfNeq: %.2f  RMS: %.2f Gap: %.2f\n",
156                pOrg->dLat, pOrg->dLon, pOrg->dZ, pOrg->dT, pOrg->dAff, pOrg->dAffGap, pOrg->dAffNumArr, pOrg->dRms, pOrg->dGap,
157                pOrg->dRms, pOrg->nEq);
158
159    for(; i<nPck; i++)
160    {
161                  pck = pPick[i];
162      CDebug::Log(DEBUG_MINOR_INFO,
163                  "%3d) %-4s %-6s %5.1f %3.0f %3.0f %5.2f %6.1f %4.1f %3.1f/%3.1f/%3.1f  %s\n",
164                                    i, pck->sSite, pck->sPhase,
165                                    pck->dDelta, pck->dAzm, pck->dToa * RAD2DEG, pck->tRes, pck->dT - pOrg->dT, 
166                  pck->dAff, pck->dAffDis, pck->dAffRes, pck->dAffPPD, pck->sTag
167                 ); 
168    }
169
170    // if the full was better, then reprocess the parameters using the FULL
171    if(dFullAff > pOrg->dAff)
172    {
173      OPCalc_ProcessOriginSimple(pOrg, pPick, nPck);
174    }
175  }
176
177  // sort the picks by resid
178  qsort(pPick, nPck, sizeof(PICK *), ComparePickiPickAsc);
179
180  return(0);
181
182}
183int OPCalc_ProcessOriginSimple(ORIGIN * pOrg, PICK ** pPick, int nPck) 
184{
185
186  double pPickDistArray[OPCALC_MAX_PICKS];
187  double pPickAzmArray[OPCALC_MAX_PICKS];
188
189  // Initialize origin vars
190        int nEq = 0;
191  int i;
192
193        double sum = 0.0;
194  double dWeight = 0.0;
195  PICK * pck;
196  int iRetCode;
197  double dGap;
198
199  PhaseType  ptPhase;
200  PhaseClass pcPhase;
201
202  // Set default RMS to absurd level.  Shouldn't this be an affinity setting.
203  pOrg->dRms = 10000.0;
204
205
206//  CDebug::Log(DEBUG_MAJOR_INFO,"Locate(I): Org(%6.3f,%8.3f,%5.1f - %.2f)\n",
207//              pOrg->dLat, pOrg->dLon, pOrg->dZ, pOrg->dT);
208
209        for(i=0; i<nPck; i++) 
210  {
211                pck = pPick[i];
212    iRetCode = OPCalc_CalculateOriginPickParams(pck,pOrg);
213    if(iRetCode == -1)
214    {
215      CDebug::Log(DEBUG_MINOR_ERROR,"ProcessOrigin(): OPCalc_CalculateOriginPickParams  failed with "
216                                    "terminal error.  See logfile.  Returning!\n");
217      return(-1);
218    }
219    else if(iRetCode == 1)
220    {
221      // could not find TravelTimeTable entry for pOrg,pck
222      CDebug::Log(DEBUG_MINOR_INFO,"ProcessOrigin(): OPCalc_CalculateOriginPickParams() could not "
223                                    "find TTT entry for Org,Pick!\n");
224      continue;
225    }
226
227    // accumulate statistics
228    sum += pck->tRes*pck->tRes*pck->ttt.dLocWeight;
229
230    dWeight += pck->ttt.dLocWeight;
231
232    // copy the delta/azimuth of the pick to proper arrays for
233    // further calculations
234    pPickDistArray[nEq] = pck->dDelta;
235
236    // DK 021105  Change gap calculation to only use P-Class phases
237    // that way if we pull in random S, PPP, and SKiKP phases,
238    // they won't affect the gap "quality" of the solution.
239    ptPhase = GetPhaseType(pck->ttt.szPhase);
240    pcPhase = GetPhaseClass(ptPhase);
241    if(( (pcPhase == PHASECLASS_P) && (pck->dDelta < OPCALC_MAX_BELIEVABLE_DISTANCE_FOR_P_PHASE )) || 
242       (ptPhase == PHASE_PKPdf)
243      )
244      pPickAzmArray[nEq] = pck->dAzm;
245    else
246      pPickAzmArray[nEq] = -1.0;
247
248
249    // increment the number of picks used to locate the origin
250    nEq++;
251   
252    // keep the arrays from being overrun
253    if(nEq > OPCALC_MAX_PICKS)
254      break;
255  }  // end for each pick in the Pick Array
256
257  // if we didn't use atleast one pick.  Give up and go home.
258  if(nEq < 1)
259          return(1);
260
261
262  // Record the number of equations used to calc the location
263  pOrg->nEq = nEq;
264
265  // Record the residual standard deviation
266  pOrg->dRms = sqrt(sum/dWeight);
267
268  /* Calculate Pick Azimuth statistics for the Origin (Gap)
269   ****************************************************/
270  // Sort the pick-azimuth array
271  qsortd(pPickAzmArray, nEq);
272
273  // If the Pick wasn't a P-class Pick, then we gave it a negative azm,
274  // so igonre all of the negative azm picks (which are all at the beginning
275  // since the array has been sorted
276        for(i=0; i<nEq-2; i++) 
277  {
278    if(pPickAzmArray[i] >= 0.0)
279      break;
280  }
281
282  CDebug::Log(DEBUG_MAJOR_INFO,
283              "ProcessOrigin(): Calculating Gap for Org %d(%5.1f, %6.1f, %3.0f) based on %d valid picks.\n ",
284              pOrg->iOrigin, pOrg->dLat, pOrg->dLon, pOrg->dZ, nEq - i);
285  // Calculate the gap for the last-to-first wrap
286        dGap = pPickAzmArray[i] + 360.0 - pPickAzmArray[nEq-1];
287        for(; i<nEq-1; i++) 
288  {
289    // for each gap between two picks, record only if largest so far
290    if(pPickAzmArray[i+1] - pPickAzmArray[i] > dGap)
291      dGap = pPickAzmArray[i+1] - pPickAzmArray[i];
292        }
293        pOrg->dGap = dGap;
294
295  /* Calculate Pick Distance statistics for the Origin
296   ****************************************************/
297  // Sort the pick-distance array
298  qsortd(pPickDistArray, nEq);
299
300  // Calculate the minimum station distance for the origin
301        pOrg->dDelMin = pPickDistArray[0];
302
303  // Calculate the maximum station distance for the origin
304        pOrg->dDelMax = pPickDistArray[nEq-1];
305
306  // Calculate the median station distance for the origin
307        pOrg->dDelMod = 0.5*(pPickDistArray[nEq/2] + pPickDistArray[(nEq-1)/2]);
308
309
310  /* Calculate the Affinity for the Origin (and Picks)
311   ****************************************************/
312  iRetCode = OPCalc_CalcOriginAffinity(pOrg, pPick, nPck);
313  if(iRetCode)
314  {
315    CDebug::Log(DEBUG_MINOR_ERROR,
316                "Error: ProcessOrigin(): OPCalc_CalcOriginAffinity() return (%d)\n",
317                iRetCode);
318    return(iRetCode);
319  }
320     
321  return(0);
322}  // End Process Location()
323
324
[2176]325int OPCalc_CalcOriginAffinity(ORIGIN * pOrg, PICK ** ppPick, int nPck) 
326{
327//      PICK *pck;
328        int i;
329  int iRetCode;
330  double dAffinitySum = 0.0;  // Sum of usable-pick affinities
331  int    iAffinityNum = 0;    // Number of usable-picks
332
333  if(nPck < 1)
334    return(-1);
335
336  // Log the calculated values
[3212]337        CDebug::Log(DEBUG_MAJOR_INFO,"Affinity(): dDelMin:%.2f dDelMod:%.2f dDelMax:%.2f\n",
[2176]338                pOrg->dDelMin, pOrg->dDelMod, pOrg->dDelMax);
339
340        // Gap affinity
341                // DK CHANGE 011504.  Changing the range of the
342    //                    dAffGap statistic
343    //         pOrg->dAffGap = 1.0+OPCalc_BellCurve(pOrg->dGap/360.0);
344    // dAffGap range is 0.0 - 2.0
[3212]345  pOrg->dAffGap = 4 * OPCalc_BellCurve(pOrg->dGap/360.0); 
346        if(pOrg->nEq < OPCALC_nMinStatPickCount && pOrg->dAffGap < 1.0)
[2176]347                pOrg->dAffGap = 1.0;
[3212]348
[2176]349//    pOrg->dAffGap = 3 * OPCalc_BellCurve(pOrg->dGap/360.0);
350
351  if(pOrg->dAffGap > 2.0)
352    pOrg->dAffGap = 2.0;
353  else if(pOrg->dAffGap < 0.5 && pOrg->dAffGap > 0.1)
354    pOrg->dAffGap = 0.5;  // don't let AffGap drop below 0.5 unless
355                          //the Gap is AWFUL!!!!
356
357
358
359        // nEq affinity (macho statistic)
360        if(pOrg->nEq < OPCALC_nMinStatPickCount)
361                pOrg->dAffNumArr = 1.0;
362        else
[3212]363                pOrg->dAffNumArr = log10((double)(pOrg->nEq));
[2176]364
365        for(i=0; i<nPck; i++) 
366  {
367    iRetCode = OPCalc_CalcPickAffinity(pOrg, ppPick[i]);
368          CDebug::Log(DEBUG_MINOR_INFO,"Affinity(): Pck[%d] Delta:%.2f Aff:%.2f\n", i, ppPick[i]->dDelta, ppPick[i]->dAff);
369
370    if(iRetCode < 0)
371    {
372          CDebug::Log(DEBUG_MINOR_ERROR,
373                  "OPCalc_CalcOriginAffinity(): Fatal Error(%d) calculating"
374                  " affinity for pick[%d]\n",
375                  iRetCode,i);
376      return(-1);
377    }
378    if(iRetCode > 0)
379      continue;
380
381    dAffinitySum+=ppPick[i]->dAff;
382    iAffinityNum++;
383  }
384
385  if(iAffinityNum > 0)
386  {
387    pOrg->dAff = dAffinitySum / iAffinityNum;
[3212]388
389    /* now lets do something more advanced. 
390       We are trying to calculate the "quality" of this origin. 
391       We're not trying to recalculate where the origin should be,
392       only trying to caculate a "quality" number.
393       Under the current basic calculations, we can run into problems
394       where a couple of random picks can have residuals close enough in
395       to associate, but far enough out to lower the quality of the origin.
396       A problem is, that as soon as you remove a pick from the origin
397       calculation, you have to do a lot of quality recalculation.  Removing
398       one pick potentially affects both AffGap affinity component as well
399       as the AffNumArr.
400       What we'd like is to find some sort of peak quality value for the Origin,
401       (a reverse saddle).  We'll try to get there by taking the full result
402       (all-picks), and remove the picks with the lowest affinity until we hit
403       a point where the origin affinity starts to decrease. 
404       Unfortunately this can't all be done in the OPCalc routines, because the
405       Gap is currently calculated in the Glock module ONLY.
406       DK 0920'06
407      ************************************************************************/
408
409    for(i=0; i<nPck; i++) 
410    {
411    }
412
413
414
[2176]415  }
416  else
417  {
418    pOrg->dAff = 0;
419  }
420  return(0);
421}  // end OPCalc_CalcOriginAffinity()
422
423
424int OPCalc_CalcPickAffinity(ORIGIN * pOrg, PICK * pPick) 
425{
426
427  double dDistRatio;  // ratio between distance of current pick and median distance
428  double dMedianDistEst;
429
430  int iRetCode;
431
432        // Do full calculations for picks not already associated with the
433  // origin.  (distance, azimuth, traveltime, TakeOffAngle)
434        if(pPick->iState == GLINT_STATE_WAIF) 
435  {
436    iRetCode = OPCalc_CalculateOriginPickParams(pPick,pOrg);
437    if(iRetCode)
438      return(iRetCode);
439        }
440
441        // Residual affinity
442  // Range   0.0 - 2.0  where 0.0 is no residual and 2.0 is
443  //         any absolute residual >= 10 seconds (OPCALC_dResidualCutOff)
444        pPick->dAffRes = 2.0 * OPCalc_BellCurve(pPick->tRes/(pPick->ttt.dResidWidth));
445
446        // Distance affinity
447  if(OPCALC_fAffinityStatistics & AFFINITY_USE_DISTANCE)
448  {
449
450    // hack in a calculation that increases the value used for the median distance if the number of
451    // phases is large.
452    if(pOrg->nPh > (pOrg->dDelMod*2.0))
453      dMedianDistEst = (double) (pOrg->nPh/2);
454    else
455      dMedianDistEst = pOrg->dDelMod;
456   
457    // Calculate the distance of this pick relative
458    // to the distance of median pick used in the solution
459    if(pOrg->nEq < OPCALC_nMinStatPickCount) 
460    {
461        dDistRatio = pPick->dDelta / (OPCALC_dCutoffMultiple*dMedianDistEst * 1.5);
462    } 
463    else 
464    {
465        dDistRatio = pPick->dDelta / (OPCALC_dCutoffMultiple*dMedianDistEst);
466    }
467    // Distance Affinity
468    //       Penalize picks (based upon distance) when
469    //       they are more than dCutoffMultiple times the
470    //       median pick distance from the Origin.
471    // Range   0.0 - 2.0  where 0.0 is beyond the Modal distance cutoff and 2.0 is
472    //         a pick at the hypocenter.
473    pPick->dAffDis = 2 * OPCalc_BellCurve(dDistRatio);
474
475    if(pOrg->nEq < 4)
476      pPick->dAffDis = 1.0;
477  }
478  else
479  {
480    pPick->dAffDis = 1.0;
481  }
482
483  if(OPCALC_fAffinityStatistics & AFFINITY_USE_PPD)
484  {
485    // Pick Probability Affinity
486    // Calculate a Pick Probability Affinity statistic.
487    // Range   0.5 - 5.0 (highest observed value is 3.0)
488    // where the value is ((log 2 (X+1)) + 1) / 2
489    // where X is the number of times that a phase was picked
490    // within the Pick's travel time entry per 10000 picks.
491    pPick->dAffPPD = pPick->ttt.dAssocStrength / 2.0;
492    if(pPick->dAffPPD > 2.0)
493      pPick->dAffPPD = 2.0;
494  }
495  else
496  {
497    pPick->dAffPPD = 1.0;
498  }
499
500        // Composite affinity statistic
501        pPick->dAff = pOrg->dAffGap * pOrg->dAffNumArr * pPick->dAffRes * pPick->dAffDis * pPick->dAffPPD;
502
503  return(0);
504}  // end CalcPickAffinity()
505
506
507void OPCalc_SetTravelTimePointer(ITravelTime * pTravTime)
508{
509  pTT = pTravTime;
510}
511
512
513int OPCalc_CalculateOriginPickParams(PICK * pPick, ORIGIN * pOrg)
514{
515  TTEntry   ttt;
516  TTEntry * pttt;
517  double    dDistKM;
518
519  if(!(pTT && pPick && pOrg))
520  {
521                CDebug::Log(DEBUG_MINOR_ERROR,"OPCalc_CalculateOriginPickParams(): Null Pointer: "
522                                  "pPick(%u) pOrg(%u) pTT(%u)\n", 
523                pPick, pOrg, pTT);
524    return(-1);
525  }
526
527  if(pPick->dLat < -90.0 || pPick->dLat > 90.0 || pPick->dLon < -180.0 || pPick->dLon > 180.0)
528    return(-1);
529
530        geo_to_km_deg(pOrg->dLat, pOrg->dLon, pPick->dLat, pPick->dLon, 
531                &dDistKM, &pPick->dDelta, &pPick->dAzm);
532        pttt = pTT->TBest(pOrg->dZ, pPick->dDelta, pPick->dT - pOrg->dT, &ttt);
533        if(!pttt) 
534  {
535    pPick->bTTTIsValid = false;
536    pPick->tRes = 9999.9;
537    pPick->ttt.szPhase = Phases[PHASE_Unknown].szName;
538
539    return(1);
540  }
541
542        pPick->dTrav = pttt->dTPhase;
543        pPick->dToa = pttt->dTOA;
544        pPick->tRes = (pPick->dT - pOrg->dT) - pttt->dTPhase;
545  strcpy(pPick->sPhase, pttt->szPhase);
546  memcpy(&pPick->ttt, pttt, sizeof(TTEntry));
547  pPick->bTTTIsValid = true;
548  CDebug::Log(DEBUG_MINOR_INFO,"Iterate(): Eq[%d] d:%.2f z:%.2f tobs:%.2f azm:%.2f ttt:%d\n", 
549              pOrg->nEq, pPick->dDelta, pOrg->dZ, pPick->dT - pOrg->dT, pPick->dAzm, pttt);
550
551  return(0);
552}  // end OPCalc_CalculateOriginPickParams()
553
554
Note: See TracBrowser for help on using the repository browser.