Develop/1.0

experiments/fullInterpolationAndStack/SConstruct

@@ -56,8 +56,8 @@ ot0 = 1.0
 dt0 = 0.004
 nt0 = 500
 om0 = 3
-dm0 = 0.025
-nm0 = 160
+dm0 = 0.5
+nm0 = 9
 
 for i in range(nm0):
 

experiments/fullInterpolationAndStack/creStack.py

@@ -160,6 +160,6 @@ def pefInterpolation(
             '''
             rcat axis=3 ${SOURCES[1:%d]} |
             transp plane=23
-            ''' % offsetGatherIndex)
+            ''' % nhi)
 
 

experiments/fullInterpolationAndStackCDS/SConstruct

@@ -0,0 +1,181 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# SConstruct  (Madagascar Script)
+#
+# Purpose: Build the interpolation process twice to improve CMP sampling
+# and CRE stacking results. 
+#
+# IMPORTANT: This SConstruct uses the non-hyperbolic CRS traveltime approximation
+# with the CDS condition applied (RN=RNIP) as a stacking curve.
+#
+# Site: http://www.dirackslounge.online
+# 
+# Version 1.0
+#
+# Programer: Rodolfo A. C. Neves (Dirack) 16/03/2020
+#
+# Email: rodolfo_profissional@hotmail.com
+#
+# License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
+
+# Madagascar package
+from rsf.proj import *
+
+# CRE recipe
+from creStack import kirchoffModeling as kimod
+from creStack import pefInterpolation as pefin
+
+# Import glob python library
+import glob
+
+# Generate Gaussian reflector model and data cube
+kimod()
+
+# PEF interpolation of the data cube
+pefin('dataCube',
+    'interpolatedDataCube',
+    nm=401,
+    dm=0.025,
+    nt=1001,
+    dt=0.004,
+    nhi=161)
+
+# Do the PEF interpolation one more time
+# to increase the CMP sampling
+pefin('interpolatedDataCube',
+    'interpolatedDataCube2',
+    nm=802,
+    dm=0.0125,
+    nt=1001,
+    dt=0.004,
+    nhi=161)
+
+# CRE stacking
+# It uses Very Fast Simulated Aneeling and non hyperbolic CRS
+# to get zero offset CRS parameters (RN, RNIP and BETA) from data cube
+v0 = 1.5
+ot0 = 1.0
+dt0 = 0.004
+nt0 = 500
+om0 = 3
+dm0 = 0.5
+nm0 = 9
+
+for i in range(nm0):
+
+	parametersCube = []
+	creGatherCube = []
+	creTimeCurveCube = []
+	creGatherCubeForM0List = []
+	creTimeCurveCubeForM0List = []
+	creGatherIndex = 0
+
+	m0 = om0 + (i * dm0)
+
+	for j in range(nt0):
+
+		t0 = ot0 + (dt0 * j)
+
+		crsParameters = 'crsParameters-m0-%g-t0-%g' % (i,j)
+		creMhCoordinates = 'creMhCoordinates-m0-%g-t0-%g' % (i,j)
+		creGather = 'creGather-m0-%g-t0-%g' % (i,j)
+		creMcoordinate = 'creMcoordinate-m0-%g-t0-%g' % (i,j)
+		creTimeCurve = 'creTimeCurve-m0-%g-t0-%g' % (i,j)
+		creGatherPlot = 'crePlot-m0-%g-t0-%g' % (i,j)
+
+		# Very Fast Simulated Aneelling Global Optimization (VFSA)
+		Flow(crsParameters,'dataCube',
+			'''
+			vfsacrsnh m0=%g v0=%g t0=%g verb=y repeat=3
+			''' % (m0,v0,t0))
+
+		# Calculate CRE trajectory
+		Flow(creMhCoordinates,['interpolatedDataCube2',crsParameters],
+			'''
+			cretrajec verb=y m0=%g param=${SOURCES[1]} |
+			put unit1="Offset" label1="Km"
+			''' % (m0))
+
+		#Get CRE Gather from interpolated Data Cube
+		Flow([creGather,creMcoordinate],['interpolatedDataCube2',creMhCoordinates],
+			'''
+			getcregather verb=y cremh=${SOURCES[1]} m=${TARGETS[1]} |
+                        put label1="Time" unit1="s" label2="Offset" unit2="km" 
+			''')
+
+		# Calculate CRE traveltime curve t(m,h)
+		Flow(creTimeCurve,[creMcoordinate, crsParameters],
+			'''
+			getcretimecurve param=${SOURCES[1]} cds=y t0=%g m0=%g v0=%g verb=y |
+                        put label1="Offset" unit1="Km" 
+			''' % (t0,m0,v0))
+
+		parametersCube.append(crsParameters)
+		creGatherCube.append(creGather)
+		creTimeCurveCube.append(creTimeCurve)
+		creGatherIndex = creGatherIndex + 1
+
+	# Concatenate cre gathers for each m0
+	creGatherCubeForM0 = "creGatherCube-m0-%i" % i
+	creGatherCubeForM0List.append(creGatherCubeForM0)
+	Flow(creGatherCubeForM0,creGatherCube,
+		'''
+		rcat axis=3 ${SOURCES[1:%d]}
+		''' % nt0)
+
+	# Concatenate cre Time curves for each m0
+	creTimeCurveCubeForM0 = "creTimeCurveCube-m0-%i" % i
+	creTimeCurveCubeForM0List.append(creTimeCurveCubeForM0)
+	Flow(creTimeCurveCubeForM0,creTimeCurveCube,
+		'''
+		rcat axis=2 ${SOURCES[1:%d]}
+		'''% nt0)
+
+	# Stack one m0 per turn
+	creTimeCurveCube = []
+	creGatherCube = []
+	creGatherTrace = "creGatherTrace-m0-%i" % i
+	creTimeCurveTrace = "creTimeCurveTrace-m0-%i" % i
+	creStackedTrace = "creStackedTrace-m0-%i" % i
+
+	# Get cre Gathers organized by t0 and m0
+	Flow(creGatherTrace,creGatherCubeForM0List,
+		'''
+		put n4=1 d4=%g o4=%g d3=%g o3=%g label3=t0 unit3=s label4=m0 unit4=Km
+		''' % (dm0,m0,dt0,ot0))
+
+	# Get all the traveltime curves organized by t0 and m0
+	Flow(creTimeCurveTrace,creTimeCurveCubeForM0List,
+		'''
+		put label3="m0" unit3="Km" label2="t0" unit2="s" d2=%g o2=%g d3=%g o3=%g
+		''' % (dt0,ot0,dm0,m0))
+
+	# CRE stacking
+	Flow(creStackedTrace,[creGatherTrace,creTimeCurveTrace],
+		'''
+		crestack timeCurves=${SOURCES[1]} verb=y |
+		put label1=t0 unit1=s label2=m0 unit2=Km --out=stdout
+		''')
+
+# Build the cre stacked section
+# throughout cre stacked traces sorting
+files = glob.glob('creStackedTrace-*.rsf')
+length = len(files)
+
+sortedFiles = []
+for i in range(length):
+	string = 'creStackedTrace-m0-%i.rsf' % i
+	sortedFiles.append(string)  
+
+Flow('stackedSection',sortedFiles,
+	'''
+	rcat axis=2 ${SOURCES[1:%d]} --out=stdout
+	''' % len(files))
+
+Flow('filtStackedSection','stackedSection',
+	'''
+	bandpass fhi=20 --out=stdout
+	''')
+
+End()

experiments/fullInterpolationAndStackCDS/creStack.py

@@ -0,0 +1,165 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# creStack  (Python)
+#
+# Purpose: Recipe to the cre model, interpolation and Stacking.
+# 
+# Important!: It should be called from a SConstruct 
+#
+# Site: http://www.dirackslounge.online
+# 
+# Version 1.0
+#
+# Programer: Rodolfo A. C. Neves (Dirack) 04/03/2020
+#
+# Email: rodolfo_profissional@hotmail.com
+#
+# License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
+
+# Madagascar package
+from rsf.proj import *
+
+def kirchoffModeling(filename='dataCube'):
+    '''
+    Modeling function of a gaussian reflector
+    '''
+
+    # Modeling: Gaussian reflector in a velocity linear model
+    # velocity increases with depth with a 0.5 velocity gradient
+    Flow('gaussianReflector',None,
+         '''
+         math d1=0.01 n1=2001 o1=-5 unit1=km label1=Offset
+         output="4-3*exp(-(x1-5)^2/9)"
+         ''')
+
+    # Velocity Model
+    Flow('velocityModel','gaussianReflector',
+         '''
+         window min1=0 max1=10 |
+         spray axis=1 n=451 d=0.01 o=0 label=Depth unit=km |
+         math output="1.5+0.5*x1+0.0*x2"
+         ''')
+
+    Flow('reflectorDip','gaussianReflector','math output="2/3*(x1-5)*input" ')
+
+    # Kirchoff Modeling
+    Flow(filename,'gaussianReflector reflectorDip',
+         '''
+         kirmod cmp=y dip=${SOURCES[1]} 
+         nh=161 dh=0.025 h0=0
+         ns=401 ds=0.025 s0=0
+         freq=10 dt=0.004 nt=1001
+         vel=1.5 gradz=0.5 gradx=0.0 verb=y |
+         put d2=0.0125 label3="CMP" unit3="Km" label2="Offset" unit2="Km" label1=Time unit1=s
+         ''')
+
+def pefInterpolation(
+    dataCube,
+    interpolated,
+    nm,
+    dm,
+    nt,
+    dt,
+    nhi=1
+    ):
+    '''
+    PEF interpolation of the data cube
+    :param dataCube: filename, Seismic data cube to interpolate
+    :param interpolated: filename, Interpolated seismic data cube
+    :param nm: integer, number of CMPs in the seismic data cube
+    :param dm: float, CMP sampling
+    :param nt: integer, number of time samples
+    :param dt: float, time sampling
+    :param nhi: integer, number of constant offsets gathers to interpolate
+    '''
+
+    # Divide CMP sampling
+    dm = dm/2
+
+    # Define mask file names using input filename
+    mask1 = dataCube+'-mask1'
+    mask = dataCube+'-mask'
+    aa = dataCube+'-aa'
+    bb = dataCube+'-bb'
+    a = dataCube+'-a'
+    b = dataCube+'-b'
+    zeroTraceGather = dataCube+'-zeroedGather'
+    mask0 = dataCube+'-mask0'
+
+
+    # Build a mask to interleave zero traces with original data traces
+    Flow(aa,None,'spike n1=%i d1=%g o1=0' %(nm,dm))
+    Flow(bb,None,'spike n1=%i d1=%g o1=0 mag=0' % (nm,dm))
+    Flow(mask1,[bb, aa],
+            '''
+            interleave axis=1 ${SOURCES[1]} |
+            dd type=int
+            ''')
+
+    Flow(a,None,'spike n1=%i d1=%g o1=0' % (nm,dm))
+    Flow(b,None,'spike n1=%i d1=%g o1=0 mag=0' % (nm,dm))
+    Flow(mask,[a, b],
+            '''
+            interleave axis=1 ${SOURCES[1]} |
+            dd type=int
+            ''')
+    Flow(zeroTraceGather,b,
+            '''
+            spray axis=2 n=%i d=%g |
+            transp |
+            put label2=Offset unit2=Km label1=Time unit1=s
+            ''' %(nt,dm))
+
+    # Data Mask with double of traces in CMP (half of CMP sampling)
+    # Keep the same Time and Offset original data sampling
+    Flow(mask0,mask,
+         '''
+         spray axis=1 n=%i d=%g
+         ''' %(nt,dt))
+
+    totalPefIterations = 100
+    totalInterpolationIterations = 20
+
+    offsetGathers = []
+    for offsetGatherIndex in range(nhi):
+
+            offsetGather = dataCube+"-offsetGather-%i" % offsetGatherIndex
+            resampledOffsetGather = dataCube+"-resampledGather-%i" % offsetGatherIndex
+            interpolatedOffsetGather = dataCube+"-interpolatedGather-%i" % offsetGatherIndex
+            pefCoeficients = dataCube+"-pefCoeficients-%i" % offsetGatherIndex
+
+            Flow(offsetGather,dataCube,
+            '''
+            window n2=1 f2=%i
+            ''' % (offsetGatherIndex))
+
+            Flow(resampledOffsetGather,[offsetGather,zeroTraceGather],
+            '''
+            interleave axis=2 ${SOURCES[1]}
+            ''')
+
+            # Calculate adaptive PEF coeficients
+            Flow(pefCoeficients,[resampledOffsetGather,mask0],
+                    '''
+                    apef jump=2 a=10,2 rect1=50 rect2=2 niter=%g verb=y
+                    maskin=${SOURCES[1]}
+                    ''' % (totalPefIterations))
+
+            # Interpolation
+            Flow(interpolatedOffsetGather, 	[resampledOffsetGather,pefCoeficients,mask0,mask1],
+                    '''
+                    miss4 exact=y filt=${SOURCES[1]} mask=${SOURCES[2]} niter=%g verb=y |
+                    put d2=%g
+                    ''' % (totalInterpolationIterations,dm))
+
+            offsetGathers.append(interpolatedOffsetGather)
+
+    # Concatenate interpolated sections
+    Flow(interpolated,offsetGathers,
+            '''
+            rcat axis=3 ${SOURCES[1:%d]} |
+            transp plane=23
+            ''' % nhi)
+
+