Seismic Data Reconstruction package

SeisRec_0.1

Author:
Mostafa Naghizadeh
Copyright (c) 2006-2007 Mostafa Naghizadeh

Inroduction

This is a C Package which contains the codes for Multi-Dimensional Seismic Reconstruction. The included reconstruction methods are:

1. Minimum Weighted Norm Interpolation (MWNI).

2. Multi-Step Auto-Regressive (MSAR) reconstruction.

3. The F-K domain interpolation (Gulunay)

This is the first version of the library and may not be as robust or well documented as it should be. Please keep track of bugs or missing/confusing instructions and report them to mostafan@ualberta.ca

References

1. Naghizadeh, M. and M. D. Sacchi (2007). Multi-step auto-regressive reconstruction of seismic records. Geophysics 72 (6), V111-V118.

2. Liu, B. and M. D. Sacchi (2004). Minimum weighted norm interpolation of seismic records. Geophysics 69 (6), 1560-1568.

3. Gulunay, N. (2003). Seismic trace interpolation in the fourier transform domain. Geophysics 68 (1), 355-369.

Installation:

1. Download the package seisrec.tar.gz

2. Unzip the package to your home directory.

>gunzip seisrec.tar.gz
>tar -xf seisrec.tar
This means you will have the directory seisrec inside your home directory.

3. You need to add a line to the file .cshrc inside your home directory. Do the following steps:

> emacs .cshrc
Now add the following line to the .cshrc file
setenv SEISREC $HOME/seisrec
Close the file and and run the following command in your command line:
> source .cshrc

4. Enter seisrec directory

> cd seisrec

5. Now you can configure the package using

./seisrec> configure -s
Configure command will creat Makefile for the libray and all the examples.

6. In order to compile and make executable file for examples, just type make in the command line:

./seisrec> make
Now package is ready and you can proceed to the example directory and run the examples.

Examples:

There are several examples provided in the example folder of the package. You should be able to see the results by just running RUNME.m file inside MATLAB environment. MATLAB is just used for making the synthetic examples and visualizing the results.

Submiting an example job

In order to have in depth understanding of using the package, we examine the parameters of msar_synthetic_2d.c located inside the example folder. All of the example have the same style of submitting the job . In order to submit a new job you can follow the style of the provided examples.

1   #include <stdio.h>
2   #include "msar.h"

3   int main(void)
   {
4    int NDIM,ia,ind,iter_cg,iter_bl;
     unsigned long *nx,*nhan,nh,js,*npf;
     float fmin,fmax,*bmin,*bmax,*rzp,fmaxp;
     msar_plan ms_plan;
     mwni_plan mw_plan;
     FILE *fid;
  
5    NDIM=2;                      // Dimension of data 
     nx=lvector(1,NDIM);       // Vector for number of samples in each direction of filter 
     rzp=vector(1,NDIM);
     bmin=vector(1,NDIM-1);
     bmax=vector(1,NDIM-1);
     nhan=lvector(1,NDIM-1);
     npf=lvector(1,NDIM-1);

 6  // Number of whole and avialable traces and time smaples
     fid=fopen("data_count.dat","r");	
     fscanf(fid,"%lu",&nx[1]);
     fscanf(fid,"%lu",&nx[2]);
     fscanf(fid,"%lu",&nh);
     fclose(fid);


7   // Hanning dimension 
     nhan[1]=3;

8   // zero-padding ratio 
     rzp[1]=0;
     rzp[2]=1.5;

9   // band-limiting values 
     fmin=0.035;
     fmax=0.07;
     bmin[1]=.45;
     bmax[1]=.55;

10  // type of MWNI 
      ind=1;

11  // Number of iterations 
     iter_cg=10;
     iter_bl=3;
  
12  // msar parameters 
     js=8;
     fmaxp=.5;
     npf[1]=5;


13  // Initiating MSAR plan 
     init_msar_plan(&mw_plan,&ms_plan, NDIM , nx, nh, rzp, fmin, fmax, 
		        bmin, bmax, ind, iter_cg, iter_bl, nhan, npf, fmaxp, js); 

14  // Giving values of avialable traces 
     fid=fopen("avail_data.dat","r");	
     for(ia=1;ia<=mw_plan.nh*mw_plan.n[1];ia++)
       {
         fscanf(fid,"%f",&mw_plan.data[ia]);
       }
     fclose(fid);

15  // Giving the location of avialable tarces 
     fid=fopen("offsets.dat","r");	
     for(ia=1;ia<=mw_plan.nh;ia++)
       {
         fscanf(fid,"%lu",&mw_plan.h[ia]);
       }
     fclose(fid);

16  // Performing MSAR reconstruction 
     perform_msar(&mw_plan, &ms_plan);
  
17  // Saving values of msar reconstructed traces 
     fid=fopen("msar_rec_data.dat","w");	
     for(ia=1;ia<=mw_plan.nr*mw_plan.n[1];ia++)
       {
         fprintf(fid,"%8.5f \n",ms_plan.rdata[ia]);
       }
     fclose(fid);

18  // Saving values of mwni reconstructed traces 
     fid=fopen("mwni_rec_data.dat","w");	
     for(ia=1;ia<=mw_plan.nr*mw_plan.n[1];ia++)
       {
         fprintf(fid,"%8.5f \n",mw_plan.datar[ia]);
       }
     fclose(fid);
  
19  // Saving Number of PFs extracted for each frequency
     fid=fopen("pfc.dat","w");	
     for(ia=1;ia<=ms_plan.ntz1/2+1;ia++)
       {
         fprintf(fid,"%lu \n",ms_plan.pfc[ia]);
       }
     fclose(fid);

20  // Finalizing MWNI paln 
      finalize_msar_plan(&mw_plan,&ms_plan);
  
21   return 0;
}
Here we explain the lines indicated by numbers in above example:

1. This line includes the C library stdio.h, which is required for standard input and output management of the data.

2. Here we include the library msar.h in our example in order to apply the MSAR reconstruction. It should be noticed that the library mwni.h is included inside msar.h and as a result it would be included in our example as well.

3. This is a standard C command to start a main C code.

4. These lines define the parameters which are going to be used inside this example. In C programming all the parameters must be defined at the top of the program before appearing inside the code. Besides the ordinary parameter formats of C such as int, float, unsigned long and FILE, we also define two other formats, msar_plan and mwni_plan, which will contain all the parameters needed for our job at the later stages.

5. In this line we define the dimension of the data. For this example we have a 2D seismic section as input so we have NDIM=2. In the following lines we allocate memory for the arrays we need to pass them as input to the algorithm. Further details of each of these parameters are given at below lines.

6. In this line, the number of samples in each direction (nx[1] (time) and nx[2](space)) and also the number of available traces (nh) are given to the code. These values are already saved in a file called “data_count.dat”.

7. This line defines the length of the Hanning window used for smoothing inside the MWNI. The reconstruction algorithm for 2D data (time-space) is applied for the elements of each single frequency. Since, in this example, we have only one space direction, the frequency slices are 1D. This means the Hanning window has to be 1D as well. For 3D data (time-space-space) the Hanning window would be a 2D window.

8. In order to have an artificial-free reconstruction in the boundaries of the data, it is important to zero-pad the sides of the data. This is because an abrupt cutting of events at the boundaries of the data causes severe artifacts in the Fourier domain. In this line user determines the ratio of zero-padding for each axis of data. Zeros are padded in both end of each axis with the indicated ratios. By setting rzp[1]=0 we do not zero-pad the time direction for this example. However for the spatial direction rzp[2]=1.5 asks to padd each end of the spatial direction with 1.5 size of the number of spatial sample (tarces). For example, if nx[2] is equal to 32 then the final size of spatial samples after zero-padding would be equal to 48+32+48=128 for this specific example. When the reconstruction is done, these zeros would be cut out from boundaries (Even though there are not zeros any more). This action would prevent any kind of artificial events appearance at the sides of the data.

9. This line determines the size of the low frequency band to be reconstructed by MWNI. We represent the frequency axis with the values in the range [0,.5]. Here we chose fmin (minmum frequency to be reconstructed using MWNI) equal to 0.035 and fmax (maximum frequency to be reconstructed using MWNI) equal to 0.07. It’s the user's duty to check the data and pick the proper values for fmin and fmax. The next step is to determine how much band-limited we want each frequency slice to be. This is done by the parameters bmin[1] and bmax[1]. Their value should be at the range (0,1). Zero means complete band-limiting (forcing all the wave-numbers to be zero) and one means no band-limiting (none of the wave-numbers are forced to be zero and all of the determined by the data). While we just define the band-limiting parameters for fmin and fmax, the band-limiting value for the frequencies in between is determined by a linear interpolation scheme. Notice that both bmin and bmax have only one element for 2D data. For the case of 3D data these parameters would have two elements each of them representing each spatial direction (for more see the 3D examples provided in example folder) and so for. Interested user are urged to try value higher than 1 for bmax and see how it effects the band-limiting for the frequencies in between.

10. In order to examine the different aspects of Fourier reconstruction, one can choose various type of applying MWNI by setting the parameter ind to 0, 1 and 2. Zero means Band-Limited Fourire Reconstruction (BLFR) without reweightening, one means MWNI (Band-limiting and weigthening) and two means Sparse Reconstruction (No band-limiting).

11. The optimization of MWNI method is done using Conjugate Gradient (CG) algorithm. The number of iterations for CG algorithm is determined by iter_cg which is set to 10 for this particular example. The CG routine has to be repeated inside an iterative routine to gain new weights for the data reconstruction purposes. The number of iterations for updating the weight function is determined by iter_bl which is fixed to be 3 for this example.

12. This line shows the parameters needed for prediction filetr estimation and MSAR method. The Maximum number of jumping steps for the MSAR algorithm is determined by js. While it is defined by user, but still it depends on the number of spatial sample and also the number of elements of prediction filter. In this case if js was bigger then the one is driven by data, it would be replaced by the new calculated js inside the code. On can also determine the maximum frequency (s)he wants to do the MSAR reconstruction by setting the value of fmaxp. Again this can not be bigger than 0.5. Also the length of prediction filters is determined by npf[1]=5 in this example. Notice that for 3D data we needed to define the length of Pf by setting another direction as npf[2].

13. In this line the above tuned parameters are fed to the subroutine inti_msar_plan. This subroutine will allocate memory and initial values for all the vectors and parameters needed for the MSAR reconstruction. For the MSAR plan we initiate one msar_plan (ms_plan) and one mwni_plan (mw_plan). These two plans will save all the parameters during application of the method and there is no need to define new parameter for each subroutine inside the code. This means by just passing the plans inside each subroutine all required information will be delivered to that specific subroutine.

14. In this line the values of available traces are given to algorithm. The available traces are already saved in the file "avail_data.dat" in a lexicographic order.

15. One also needs to give the locations of the available traces in the spatial direction. The locations of traces are also is saved a file (offsets.dat) in a lexicographic order (for more than one spatial dimension).

16. This line is the main part of applying the MSAR reconstruction. The subroutine called "perform_msar" will take the values of parameters initiated inside mw_paln and ms_plan and perform the MSAR reconstruction. The desired results are saved in the proper vectors located inside these plans.

17. In this line gives the result of MSAR reconstruction are saved in the file called "msar_rec_data.dat".

18. This line saves the result of MWNI part of data in the file "mwni_rec_data.dat".

19. In this line the number of extracted prediction filters for each frequency is saved in the file "pfc.dat".

20. Like all computationally intensive C programs, one should free all the vectors which have been allocated. In this line all the allocated memories for vectors inside the MSAR and MWNI plans are freed.

21. Who cares what this line does?

After running each example you can also make postscript plot and x plot presentation of the results by running (if that specific exampl has a file called PSplot.sh in its folder and you have Seismic Unix (SU) installed in your PC):

seisrec/your_example> PSplot.sh

In order to clean the results and end up with the source files that were in the folder before you run the example do the followings in your command line of x window:

seisrec/your_example> make clean
seisrec/your_example> make
You can also run these linux commands inside your matlab window by just adding "!" before the comannds.

Bugs:

Please report any problems to mostafan@ualberta.ca
SeisRec_0.1 (c) 2006-2007 Mostafa Naghizadeh