从CSD中求二维空间谱的正确方法

Question

​

​

我试着根据上面的公式实现空间谱(见附件)。

其中kX，kY是k空间中的网格点，I‘the和j’the传感器之间的C(w，r) -交叉谱密度(这里是一个大小为ns * ns >no的矩阵。(传感器)。x，y是传感器之间的距离。( kx，ky的网格密度)

我寻找上述公式的适当python实现。我有34个传感器，它产生大小为[row*column]=[n*34]的数据。首先，我在每个传感器的数据中找到了交叉光谱密度(CSD)。然后对CSD值进行二维DFT，得到空间谱。

*)我不知道这个程序是否正确。**) python实现过程正确吗？*)另外，如果有人提供了一些相关的教程/链接，这也将对我有帮助。

​

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
import cmath

# Finding cross spectral density (CSD)
fs=500
def csdMat(data):
    rows, cols = data.shape
    total_csd = []
  
    for i in range(cols):
 
        for j in range(cols):
            f, Pxy = signal.csd(data[:,i], data[:,j], fs, nperseg=512)
            abs_csd = np.abs(Pxy)
            total_csd.append(abs_csd)                     # output as list
            csd_mat = np.array(total_csd)
    return csd_mat

## Spatial Spectra:- DFT of the csd along two dimension

def DFT2D(data):
    #data = np.asarray(data)
    dft2d = np.zeros((M,N), dtype=complex)
    for k in range(len(kx)):
        for l in range(len(ky)):
            sum_matrix = 0.0
            for m in range(M):
                for n in range(N):
                    e = cmath.exp(- 1j * ((kx[k] * dx[m]) / len(dx) + (ky[l] * dy[n]) / len(dy)))
                    sum_matrix +=  data[m,n] * e
            dft2d[k,l] = sum_matrix
    return dft2d

raw_data=np.reshape(np.random.rand(10000*34),(10000,34))

# Call the seismic array
#** Open .NPY files as an array 
#with open('res_array_1000f_131310.npy', 'rb') as f:
#    arr= np.load(f)
#raw_data = arr[0:10000, :]

#CSD of the seismic data
csd = csdMat(raw_data)
print('Shape of CSD data', csd.shape)

# CSD data of a specific frequency
csd_dat=csd[:, 11]  
fcsd = np.reshape(csd_dat, (-1, 34))
fcsd.shape

n = 34
f = 10  # frequency in Hz
c = 50  # wave speed 50, 80, 100, 200  m/s
k = 2.0*np.pi*f/c  # wavenumber
nx = n  # grid density
ny = n
kx = np.linspace(-k,k,nx)  # space vector
ky=  np.linspace(-k,k,ny)   # space vector

# Distance[Meter] between sensors 
x = [2.1,2.1,-0.7,-2.1,-2.1,-0.7,-0.7,0.6,-5.7,-8.5,-11.4,-7.7,-6.3,-3.5,-2.1,-3.4,5.4,-5.2,-8.9,-10,-10,5.4,5.4,-0.8,-3.6,-6.2,-6.8,-12.2,-17.1,-19,-18.6,-13.5,14.8,14.8]
y = [6.65,4.15,3.65,5.05,7.25,8.95,11.85,8.95,-2,-0.6,-0.9,1.25,2.9,0.9,-0.1,-1.4,9.2,5.2,4.8,6.1,8.9,13.3,17.1,17.9,13.8,-9.3,-5.2,-3.6,-3.6,-0.9,3.7,3.7,-1.8,5.7]

dx = np.array(x);  M = len(dx)
dy = np.array(y) ; N = len(dy)
X,Y = np.meshgrid(kx, ky)

dft = DFT2D(fcsd)  # Data or cross-correlation matrix
spec = dft.real    # Spectrum or 2D_DFT of data[real part]

spec = spec/spec.max()

plt.figure()
c = plt.imshow(spec, cmap ='seismic', vmin = spec.min(), vmax = spec.max(),
                 extent =[kx.min(), kx.max(), ky.min(), ky.max()],
                interpolation ='nearest', origin ='lower')
plt.colorbar(c)
plt.rcParams.update({'font.size': 18})
plt.xlabel("Wavenumber, $K_x$ [rad/m]", fontsize=18)
plt.ylabel("Wavenumber,$K_y$ [rad/m]", fontsize=18)
plt.title(f'Spatial Spectrum @10Hz', weight="bold")


#c = Wave Speed; 50, 80,100,200
cc = 2*np.pi*f /c *np.cos(np.linspace(0, 2*np.pi, 34)) 
cs = 2*np.pi*f /c *np.sin(np.linspace(0, 2*np.pi, 34))
plt.plot(cc,cs)

我想生成如下图01所示的数字

但是，通过使用改进的代码，我得到了分辨率较高的图形，如图02所示，这与图01不同。

​

我增加了另外两个数字来与图01进行比较。当考虑范围-k，k时，图如图03所示

它类似于w.r.t。XY-轴到图01，我认为这个数字是可以的，除了一些K空间的遗漏。我希望这里存在一个需要解决的问题。

在图04中，我们考虑k空间范围-20k，20k，它看起来很好，但在图01中没有相似的轴。

​

我将更新的数字如下：

有人能帮我生成图01或类似的类型吗？我对图02感到困惑。有人能帮我弄明白吗？提前谢谢。

Answer 1

在我看来你好像在放大中央叶。这也解释了为什么标度不能从0降到1。

如果我改变了这几行：

kx = np.linspace(-20*k,20*k,nx)  # space vector
ky=  np.linspace(-20*k,20*k,ny)   # space vector

然后我得到

​

​

看起来更接近你要找的东西。

为了提高分辨率，我做了一些重写，以获得这张新照片。请参阅下面更新的代码。

注意:我仍然不确定这样做是否正确。

​

​

我用的代码

# Code from https://stackoverflow.com/questions/70768384/right-method-for-finding-2-d-spatial-spectrum-from-cross-spectral-densities

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
import cmath

# Set up data
# Distance[Meter] between sensors 
x = [2.1,2.1,-0.7,-2.1,-2.1,-0.7,-0.7,0.6,-5.7,-8.5,-11.4,-7.7,-6.3,-3.5,-2.1,-3.4,5.4,-5.2,-8.9,-10,-10,5.4,5.4,-0.8,-3.6,-6.2,-6.8,-12.2,-17.1,-19,-18.6,-13.5,14.8,14.8]
y = [6.65,4.15,3.65,5.05,7.25,8.95,11.85,8.95,-2,-0.6,-0.9,1.25,2.9,0.9,-0.1,-1.4,9.2,5.2,4.8,6.1,8.9,13.3,17.1,17.9,13.8,-9.3,-5.2,-3.6,-3.6,-0.9,3.7,3.7,-1.8,5.7]

if (len(x) != len(y)):
    raise Exception('X and Y lengthd differ')

n = len(x)
dx = np.array(x);  M = len(dx)
dy = np.array(y) ; N = len(dy)

np.random.seed(12345)
raw_data=np.reshape(np.random.rand(10000*n),(10000,n))

f = 10  # frequency in Hz
c = 50  # wave speed 50, 80, 100, 200  m/s
k = 2.0*np.pi*f/c  # wavenumber
kx = np.linspace(-20*k,20*k,n*10)  # space vector
ky=  np.linspace(-20*k,20*k,n*10)   # space vector


# Finding cross spectral density (CSD)
fs=500
def csdMat(data):
    rows, cols = data.shape
    total_csd = []
  
    for i in range(cols):
        for j in range(cols):
            f, Pxy = signal.csd(data[:,i], data[:,j], fs, nperseg=512)
            #real_csd = np.real(Pxy)
            total_csd.append(Pxy)                     # output as list
            
    return np.array(total_csd)

## Spatial Spectra:- DFT of the csd along two dimension

def DFT2D(data):
    #data = np.asarray(data)
    dft2d = np.zeros((len(kx),len(ky)), dtype=complex)
    for k in range(len(kx)):
        for l in range(len(ky)):
            sum_matrix = 0.0
            for m in range(M):
                for n in range(N):
                    e = cmath.exp(- 1j * ((kx[k] * dx[m]) / len(dx) + (ky[l] * dy[n]) / len(dy)))
                    sum_matrix +=  data[m,n] * e
            dft2d[k,l] = sum_matrix
    return dft2d


# Call the seismic array
#** Open .NPY files as an array 
#with open('res_array_1000f_131310.npy', 'rb') as f:
#    arr= np.load(f)
#raw_data = arr[0:10000, :]

#CSD of the seismic data
csd = csdMat(raw_data)
print('Shape of CSD data', csd.shape)

# CSD data of a specific frequency
csd_dat=csd[:, 11]  
fcsd = np.reshape(csd_dat, (-1, n))

dft = DFT2D(fcsd)  # Data or cross-correlation matrix
spec = np.abs(dft) #dft.real    # Spectrum or 2D_DFT of data[real part]

spec = spec/spec.max()

plt.figure()
c = plt.imshow(spec, cmap ='seismic', vmin = spec.min(), vmax = spec.max(),
                 extent =[kx.min(), kx.max(), ky.min(), ky.max()],
                interpolation ='nearest', origin ='lower')
plt.colorbar(c)
plt.rcParams.update({'font.size': 18})
plt.xlabel("Wavenumber, $K_x$ [rad/m]", fontsize=18)
plt.ylabel("Wavenumber,$K_y$ [rad/m]", fontsize=18)
plt.title(f'Spatial Spectrum @10Hz', weight="bold")

Answer 2

根据上述方程，空间谱的脚本是较好的匹配如下。本文对"DFT2D()“函数进行了修正，使其满足方程的要求。

    def DFT2D(data):

        P=len(kx)
        Q=len(ky)
        dft2d = np.zeros((P,Q), dtype=complex)
        for k in range(P):
            for l in range(Q):
                sum_matrix = 0.0
                for m in range(M):
                    for n in range(N): #
                        e = cmath.exp(-1j*(float(kx[k]*(dx[m]-dx[n])+ float(ky[l]*(dy[m]-dy[n])))))#* cmath.exp(-1j*w*t[n]))
                        sum_matrix += data[m, n] * e
                        #print('sum matrix would be', sum_matrix)
                #print('sum matrix would be', sum_matrix)
                dft2d[k,l] = sum_matrix
         return dft2d