 Timestamp:
 12/27/07 13:05:23 (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/seismic_processing/glass/src/modules/OPCalc/opcalc.cpp
r2176 r3212 8 8 #include <opcalc.h> 9 9 #include <opcalc_const.h> 10 11 extern "C" { 12 void __cdecl qsortd (double *base, unsigned num); 13 } 14 15 16 int OPCalc_ProcessOriginSimple(ORIGIN * pOrg, PICK ** pPick, int nPck); 17 10 18 11 19 // Lib Global variables  File Scope … … 36 44 37 45 46 int __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 61 int __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 76 int __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 92 int 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], nPcki); 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 } 183 int 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 PClass 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 pickazimuth array 271 qsortd(pPickAzmArray, nEq); 272 273 // If the Pick wasn't a Pclass 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<nEq2; 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 lasttofirst wrap 286 dGap = pPickAzmArray[i] + 360.0  pPickAzmArray[nEq1]; 287 for(; i<nEq1; 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 pickdistance 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[nEq1]; 305 306 // Calculate the median station distance for the origin 307 pOrg>dDelMod = 0.5*(pPickDistArray[nEq/2] + pPickDistArray[(nEq1)/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 38 325 int OPCalc_CalcOriginAffinity(ORIGIN * pOrg, PICK ** ppPick, int nPck) 39 326 { … … 48 335 49 336 // Log the calculated values 50 CDebug::Log(DEBUG_M INOR_INFO,"Affinity(): dDelMin:%.2f dDelMod:%.2f dDelMax:%.2f\n",337 CDebug::Log(DEBUG_MAJOR_INFO,"Affinity(): dDelMin:%.2f dDelMod:%.2f dDelMax:%.2f\n", 51 338 pOrg>dDelMin, pOrg>dDelMod, pOrg>dDelMax); 52 339 … … 56 343 // pOrg>dAffGap = 1.0+OPCalc_BellCurve(pOrg>dGap/360.0); 57 344 // dAffGap range is 0.0  2.0 58 if(pOrg>nEq < OPCALC_nMinStatPickCount) 345 pOrg>dAffGap = 4 * OPCalc_BellCurve(pOrg>dGap/360.0); 346 if(pOrg>nEq < OPCALC_nMinStatPickCount && pOrg>dAffGap < 1.0) 59 347 pOrg>dAffGap = 1.0; 60 else 61 pOrg>dAffGap = 4 * OPCalc_BellCurve(pOrg>dGap/360.0); 348 62 349 // pOrg>dAffGap = 3 * OPCalc_BellCurve(pOrg>dGap/360.0); 63 350 … … 74 361 pOrg>dAffNumArr = 1.0; 75 362 else 76 pOrg>dAffNumArr = log10(( float)(pOrg>nEq));363 pOrg>dAffNumArr = log10((double)(pOrg>nEq)); 77 364 78 365 for(i=0; i<nPck; i++) … … 99 386 { 100 387 pOrg>dAff = dAffinitySum / iAffinityNum; 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 (allpicks), 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 101 415 } 102 416 else
Note: See TracChangeset
for help on using the changeset viewer.