suresstat.c

su 的源代码库

			for (j=1; j<=cmpntr[i]; j++)
				for (it=1; it<=nt; it++)
					model[i][it] += data[cmp_loc[i][j]][it];

/* normalize - no, I cannot normalize the individual traces in the supertrace */
			/* j=(1>cmpntr[i]) ? 1 : cmpntr[i];	*/
			/* for (it=1; it<=nt; it++)		*/
			/*	model[i][it] = model[i][it]/j;	*/
		}


		/* Loop over shots */
		for (ishot=1; ishot<=ns; ishot++){
		fprintf(stderr,"\nstarting shot loop, shot=%i\n", ishot);
		if (shotntr[ishot] > 1) {

			/* construct shot and cmp super traces */
			fprintf(stderr,"creating supertraces\n");
			for (it=1; it<=ntotal_super; it++) g_trace[it] = 0.;
			for (it=1; it<=ntotal_super; it++) model_trace[it] = 0.;

			for (ichan=1; ichan<=shotntr[ishot]; ichan++){
				j = (ichan-1)*nt_super + ntpick;
				l = header[shot_loc[ishot][ichan]][1];
				/* j=location in supertrace; l=cdp number */
					for (it=1; it<=nt; it++){
						g_trace[j+it] = data[shot_loc[ishot][ichan]][it];
						model_trace[j+it] = model[l][it];
				}
                        }

			/* subtract shot supertrace from cmp supertrace if desired */
			if (sub == 1) fprintf(stderr,"subtracting supertraces\n");
			if (sub == 1) for (it=1; it<=ntotal_super; it++) 
                                        model_trace[it] = model_trace[it] - g_trace[it];


			/* cross-correlate super trace 1 with super trace 2 */
			fprintf(stderr,"crosscorrelating supertraces\n");

			for (i=1; i<=ntout; i++){
				corr_trace[i] = 0.0;
				tlag=i-ntpick-1;
				for (j=ntpick; j<=(ntotal_super-ntpick); j++){
					corr_trace[i] = corr_trace[i] + g_trace[j+tlag]*model_trace[j];
				}
			}

			/* for (i=1; i<=ntout; i++) fprintf(stderr,"i=%i, corr_trace(i)=%f\n",i,corr_trace[i]); */



			/* pick cross-correlation peak */ 
			/*window(corr_trace, ntout, ntpick, pick_tr);*/
			tlag = max(corr_trace, mode, perc, ntout);

			/* remember initial estimation and total correction */
			sstat[ishot] = tlag;
			tsstat[ishot] += tlag;
			fprintf(stderr,"finished shot %i, sstat=%i, tsstat=%f (samples)\n", ishot,sstat[ishot],tsstat[ishot]);

		}
		}

		/* end shot statics loop */

		/* correct traces for shot statics (use g_trace as temp trace) */ 
		for (i=1; i<=ntraces; i++){
			ishot=header[i][3];
			for (it=1; it<=nt; it++) g_trace[it] = 0.0;
			if(sstat[ishot] > 0) {
				for (it=1; it<=(nt-sstat[nt]+1); it++) g_trace[it] = data[i][it+sstat[ishot]];
				}
			if(sstat[ishot] < 0) {
				for (it=(sstat[ishot]+1); it<=nt; it++) g_trace[it-sstat[ishot]] = data[i][it];
				}
			for (it=1; it<=nt; it++) data[i][it] = g_trace[it];
			}

		/* construct CMP stack of corrected traces */
		fprintf(stderr,"constructing new CMP stack using shot-corrected traces\n");
		for (i=1; i<=nc; i++){

			for (it=1; it<=nt; it++) model[i][it] = 0.;

			for (j=1; j<=cmpntr[i]; j++)
				for (it=1; it<=nt; it++)
					model[i][it] += data[cmp_loc[i][j]][it];

	    		/* normalize - no, I cannot normalize the individual traces in the supertrace */
			/* j=(1>cmpntr[i]) ? 1 : cmpntr[i];	*/
			/* for (it=1; it<=nt; it++)		*/
			/*	model[i][it] = model[i][it]/j;	*/
		}



		/* Loop over receivers */
		fprintf(stderr,"\n\nstarting reciever loop\n");
		for (irec=1; irec<=nr; irec++){
                if (recntr[irec] > 1) {

                        /* construct receiver and cmp super traces */
                        fprintf(stderr,"creating supertraces\n");
                        for (it=1; it<=ntotal_super; it++) g_trace[it] = 0.;
                        for (it=1; it<=ntotal_super; it++) model_trace[it] = 0.;

                        for (ichan=1; ichan<=recntr[irec]; ichan++){
                                /* j=location in supertrace; l=cdp number */
                                j = (ichan-1)*nt_super + ntpick;
                                l = header[rec_loc[irec][ichan]][1];
                                        for (it=1; it<=nt; it++){
                                                g_trace[j+it] = data[rec_loc[irec][ichan]][it];
                                                model_trace[j+it] = model[l][it];
                                }
                        }
                        /* for (i=1; i<=nt; i++) fprintf(stderr,"i=%i, model_trace(i)=%f\n",i,model_trace[i]); */
                        /* for (i=1; i<=nt; i++) fprintf(stderr,"i=%i, g_trace(i)=%f\n",i,g_trace[i]); */


                        /* subtract receiver supertrace from cmp supertrace if desired */
                        if (sub == 1) fprintf(stderr,"subtracting supertraces\n");
                        if (sub == 1) for (it=1; it<=ntotal_super; it++) 
                                        model_trace[it] = model_trace[it] - g_trace[it];


                        /* cross-correlate super trace 1 with super trace 2 */
                        fprintf(stderr,"crosscorrelating supertraces\n");

                        for (i=1; i<=ntout; i++){
                                corr_trace[i] = 0.0;
                                tlag=i-ntpick-1;
                                for (j=ntpick; j<=(ntotal_super-ntpick); j++)
					corr_trace[i] = corr_trace[i] + g_trace[j+tlag]*model_trace[j];
                        }

                        /* for (i=1; i<=ntout; i++) fprintf(stderr,"i=%i, corr_trace(i)=%f\n",i,corr_trace[i]); */



                        /* pick cross-correlation peak */
                        /*window(corr_trace, ntout, ntpick, pick_tr);*/
                        tlag = max(corr_trace, mode, perc, ntout);

                        /* remember initial estimation and total correction */
                        rstat[irec] = tlag;
                        trstat[irec] += tlag;
                        fprintf(stderr,"finished receiver %i, rstat=%i, trstat=%f (samples)\n", irec,rstat[irec],trstat[irec]);

	fprintf(stderr,"finished one receiver\n");

              }
                }

                /* end receiver statics loop */

                /* correct traces for receiver statics (use g_trace as temp trace) */
                for (i=1; i<=ntraces; i++){
                        irec=header[i][2];
                        for (it=1; it<=nt; it++) g_trace[it] = 0.0;
                        if(rstat[irec] > 0) {
                                for (it=1; it<=(nt-rstat[nt]+1); it++) g_trace[it] = data[i][it+rstat[irec]];
                                }
                        if(rstat[irec] < 0) {
                                for (it=(rstat[irec]+1); it<=nt; it++) g_trace[it-rstat[irec]] = data[i][it];
                                }
                        for (it=1; it<=nt; it++) data[i][it] = g_trace[it];
                        }

                /* construct CMP stack of corrected traces */
		if(iter<niter){
                fprintf(stderr,"constructing new CMP stack using shot-corrected traces\n");
                for (i=1; i<=nc; i++){

                        for (it=1; it<=nt; it++) model[i][it] = 0.;

                        for (j=1; j<=cmpntr[i]; j++)
                                for (it=1; it<=nt; it++)
                                        model[i][it] += data[cmp_loc[i][j]][it];

                        /* normalize - no, I cannot normalize the individual traces in the supertrace */
                        /* j=(1>cmpntr[i]) ? 1 : cmpntr[i];     */
                        /* for (it=1; it<=nt; it++)             */
                        /*      model[i][it] = model[i][it]/j;  */
		}
                }
	}
	/* end iterations */

	/* output final statics */ 
	for (it=1; it<=ns; it++) tsstat[it] = tsstat[it]*mdt;
        fps    = efopen(ssol,"wb");
	efwrite(tsstat,sizeof(float),ns,fps);
        efclose(fps);

	for (it=1; it<=nr; it++) tsstat[it] = tsstat[it]*mdt;
        fpr    = efopen(rsol,"wb");
	efwrite(trstat,sizeof(float),nr,fpr);
        efclose(fpr);

	return EXIT_SUCCESS;
}

void window (float *trace, int nt, int nnt, float *ntrace)
{
	int j;	
	int ft;	
	
	if (nnt%2 == 0) ++nnt;
	ft = (nt - 1)/2 - (nnt - 1)/2;
	for (j = 0; j < nnt; j++) {
		ntrace[j] = trace[j+ft];
	}
}

int max (float *trace, int mode, float perc, int nt)
{
	float maxamp = 0.;	/* max. amplitude sample rate */
	int globmax = 0;	/* sample of global max value */
	int locmax = 0;		/* sample of local max value */
	int zero; 		/* zero-lag sample */
	int lag = 0; 		/* lag picked */
	int j;			/* counter */

	/* Set parameters */
	if (nt%2 == 0) ++nt;
	zero = (nt+1)/2;
	globmax = (nt+1)/2;
	perc = 1.0-(perc/100);

	/* determine global max for each trace*/
	for (j = 1; j <= nt; j++) {
		if (trace[j] > maxamp) {
			maxamp = trace[j];
			globmax = j;
		}
	}
	if (mode == 1) {
	/* determine max within perc% of global max, but closer to zero lag*/
		for (j = 1; j <= nt; j++) {
			if ((trace[j] > perc*maxamp) && 
			  (ABS(j-zero) < ABS(locmax-zero))) 
				locmax = j;
		}
	}

	/* Output the result */
	switch(mode) {
		case 0:
			lag = globmax - zero;
			break;
		case 1:
			lag = locmax - zero;
			break;
	}
	return lag;
}