RosenbergerAlgorithm.py 文件源码

def rosenberger(dataX, dataY, dataZ, lambda_):
    """
    Separate P and non-P wavefield from 3-component data.

    Return a two set of 3-component traces.
    """
    # Construct the data matrix
    A = np.vstack((dataZ, dataX, dataY))

    # SVD of the first 3 samples:
    U, D, V = np.linalg.svd(A[:, 0:3])

    # Get rid of the smallest eigenvalue
    D = D[0:2]
    D = np.diagflat(D)
    U = U[:, 0:2]

    save_U = np.zeros(len(dataX))
    save_U[0] = abs(U[0, 0])

    Dp = np.zeros((3, len(dataX)))
    Ds = np.zeros((3, len(dataX)))
    Dp[:, 0] = abs(U[0, 0]) * A[:, 2]
    Ds[:, 0] = (1 - abs(U[0, 0])) * A[:, 2]

    # Loop over all the values
    for i in range(1, A.shape[1]):
        d = A[:, i]
        U, D = _update_(U, D, d, lambda_)

        Dp[:, i] = abs(U[0, 0]) * d
        Ds[:, i] = (1-abs(U[0, 0])) * d

        save_U[i] = abs(U[0, 0])

    return Dp, Ds, save_U