dcprocessing.py 文件源码

def phase_shift( shotGather, num_vel=2048, min_frequency=5, max_frequency=100, min_velocity=1, max_velocity=1000 ):

    # Ensure that min_velocity is greater than zero for numerical stability
    if min_velocity < 1:
        min_velocity = 1

    # Frequency vector
    freq = np.arange(0, shotGather.fnyq, shotGather.df)

    # FFT of timeHistories (Equation 1 of Park et al. 1998).....................
    U = np.fft.fft(shotGather.timeHistories, axis=0)

    # Remove frequencies above/below specificied max/min frequencies and downsample (if required by zero padding)
    fminID = np.argmin( np.absolute(freq-min_frequency) )
    fmaxID = np.argmin( np.absolute(freq-max_frequency) )
    freq_id = range(fminID,(fmaxID+1), shotGather.multiple)
    freq = freq[freq_id]
    U = U[freq_id, :]

    # Trial velocities
    v_vals = np.linspace( min_velocity, max_velocity, num_vel )

    # Initialize variables
    v_peak = np.zeros( np.shape(freq) )
    V = np.zeros( (np.shape(v_vals)[0], len(freq)) )
    pnorm = np.zeros( np.shape(V) )

    # Transformation ...........................................................
    # Loop through frequencies
    for c in range( len(freq) ):
        # Loop through trial velocities at each frequency
        for r in range( np.shape(v_vals)[0] ):
            # Set power equal to NaN at wavenumbers > kres
            if v_vals[r] < (2*np.pi*freq[c]/shotGather.kres):
                V[r,c] = float( 'nan' )
            # (Equation 4 in Park et al. 1998)
            else:
                V[r,c] = np.abs( np.sum( U[c,:]/np.abs(U[c,:]) * np.exp( 1j*2*np.pi*freq[c]*shotGather.position / v_vals[r] ) ) )

        # Identify index associated with peak power at current frequency
        max_id = np.nanargmax( V[:,c] )
        pnorm[:,c] = V[:,c] / V[max_id,c]
        v_peak[c] = v_vals[max_id]

    # Create instance of DispersionPower class..................................
    dispersionPower = dctypes.DispersionPower( freq, v_peak, v_vals, 'velocity', shotGather.kres, pnorm )
    return dispersionPower 



#*******************************************************************************
# Slant-stack transform