地震波场的MATLAB解析解

Question

我一直在试着求波在a中的解析解均匀，无限媒体。对于给定的源$S(t)$，波场可以在给定速度$v$的距离$r$处计算

$$W(r)=F^{-1} [-i \pi S(\ ω)H_0^{(2)} (kr)]$$

这里$F^{-1}$表示傅里叶反变换。

$S(\omega)$是$S(t)$的频域表示

$H_0^{(2)} (kr)$是零阶第二类Hankel函数。

k是波数$(=\frac{\omega}{v})$。

我跟着这个链接用Python演示代码，但我无法在MATLAB中实现相同的代码。我也怀疑它是如何运作的。由于傅里叶变换有正频率和负频率(即$-f_{max}/2

我写了下面的MATLAB代码，但得到的波场所有条目都是NaN..!

Edit1

作为一个小修改，我强制Hankel函数的第一项为零(因为它在实部包含一个NaN)。现在它工作得很好，但它产生了两个峰。我认为这类似于频移效应，正如FFT后观察到的正频率和负频率，但它不是。在这种情况下，振幅是相反的。有人能解释一下吗?

Edit2

如果我假设解决方案中只有第一个峰值是正确的，那么第二个峰值的出现是因为某种影响(我不知道)。问题来了，当我改变dt(例如$dt \in 10^{-4}[1,2,3,4，....，10]$)解的形状和振幅会发生变化。

任何帮助都是感激的。

clc;清除所有;关闭所有;vel=2000;介质速度% r=1000;波场观测距离%最终小波出现时间，t=r/vel=0.5秒t= 1;%总时间dt=.0001;%时间步长f0=25;N = round(T/dt);%样本数目t0 = 5/(√(2)*pi*f0); % Zero time shift t = dt*(0:N-1); % time vector tau = t-t0; % time vector shifted by t0 src = (1 -2* tau.*tau * f0^2 * pi^2).*exp(-tau.^2 * pi^2 * f0^2); %source wavelet %%FT of signal nfft=2^nextpow2(N); % no of fft points fmax=1/dt; % max frequency freq= fmax*(-nfft/2:nfft/2-1)/nfft; % Freq vector: -Fmax/2 < f < +Fmax/2 sw=fft(src,nfft)/nfft; % fourier transform of wavelet, S(w) sw_shift=fftshift(sw); % shifting to make S(w) amplitude corrosponding to above "freq" vector amp= abs(sw_shift); % amplitude of the S(w) %% Green Function %(I have taken only positive freq. since hankel func defined only for positive numbers. % Suggestions please..May be I am wrong..!) freq_new= fmax*(0:nfft/2-1)/nfft; w=2*pi*freq_new; const = -1i*pi; H02 = besselh(0,2,w*r/vel); H02(1)=eps; GF = const*H02; %% Calculating Wavefield %I have made Green Function (GF) symmetric and then multipled with source spectrum. GF_new=[fliplr(GF),GF]; WF= ifft(GF_new.*sw); t_new=dt*(0:length(WF)-1); %% Plot figure(); % Wavelet subplot(2,2,1); plot(t,src); title('\bf{Source wavelet}'); xlabel('Time(sec)'); ylabel('Amplitude ') % Frequencies subplot(2,2,2); plot(freq,amp); title('\bf{Source freq spec}'); xlabel('Frequencies(Hz)'); ylabel('Amplitude ') % Gneens Fkt subplot(2,2,3); plot(freq_new,abs(GF)); title('\bf{Green Function}'); xlabel('Frequencies(Hz)'); ylabel('Amplitude ') % Wafefield subplot(2,2,4); plot(t_new,real(WF)); title('\bf{Trace}'); xlabel('Time(sec)'); ylabel('Amplitude')

Answer 1

我知道很久以前有人问过这个问题。但本着Stack Exchange的精神，江南电子竞技平台我可能会把我的答案发给未来的用户。

格林函数(即解析解)

首先要注意的是波动方程在(1)源类型和(2)维数变化时的各种解。

源类型

假设声波在齐次域，这意味着你要么在解:$ $ \压裂{\部分^ 2 p (x, t)}{\部分t ^ 2} = c ^ 2 (x) \离开(\微分算符^ 2 p (x, t) + f (x, t) \右),{1}$ $ \标签($ p $是压力,美元加元纵波速度和$ f $源函数)，或者你正在求解:$ $ \左\{{数组}{rl} \ \开始压裂{\部分p (x, t)}{\部分t} & = - c ^ 2 (x) \境(\微分算符\ cdot v (x, t) + g (x, t) \境),\ \ \压裂{\部分v (x, t)}{\部分t} & = - \微分算符p (x, t)。数组{}\ \端。{2} $ $ \标签(1)与(2)之间的差异是时间差异:$$ f(x,t) = \frac{\partial}{\partial t}g(x,t)。$ $这是源函数的一个重要区别!

格林函数方程

将式(1)重申为以下形式是标准的:$ $ \离开(\压裂{1}{c ^ 2 (x)} \压裂{\部分^ 2}{\部分t ^ 2} - \微分算符^ 2 \右)p (x, t) = f (x, t)。{3} $ $ \标签把这个方程变换到频域$\partial^2/\partial t^2 \到-\²$，以及在哪里\ω^ 2美元/ c ^ 2 = k ^ 2美元我们获得:$ $ \离开(- k ^ 2 - \微分算符^ 2 \右)p (x) \ω)= f (x) \ω),{3 b} $ $ \标签或$ $ \离开(k ^ 2 + \微分算符^ 2 \右)p (x) \ω)= - f (x) \ω),{3 c} $ $ \标签我们可以咨询，例如，https://www3.nd.edu/~atassi/Teaching/AME%2060633/Notes/greens.pdf求格林函数(其中$ f (x) \ω)= \δ(x-x_s)美元与间的美元源位置)。

格林的功能

我复制一下结果。对于(1)型方程，我们有格林函数G美元：$ ${数组}{lrl} \ \开始文本{一维:}& G & =−\压裂{我}{2 k} e ^{−翼\ lvert x−间\ rvert}, {2 - d:} \ \ \文本& G & =−\压裂{我}{4}H_0 ^ {(2)} (k \ lvert x−间\ rvert),文本{3 d:} \ \ \ & G & = \压裂{1}{4π\}\压裂{e ^{−翼\ lvert x−间\ rvert}} {\ lvert x−间\ rvert}。\{数组}$ $的确,美元H_0 ^ {(2)} $是汉克尔函数，在MATLAB中我们可以简单地写出来besselh (0, 2 k * r)．

MATLAB实现

我对你的代码做了以下修改:

1. 频率向量可以从在(0,\ $ f \压裂{1}{\δt})美元而不是你是怎么做到的₂\美元(- \压裂{1}{2 \δt} \压裂{1}{2 \δt})美元．这样，汉克尔函数只接收正频率。此外，这正是FFT产生输出的顺序。我于是也去掉了fftshift操作。
2. 我把格林函数换成了二维情况，也就是说，我改变了const我*π= 1来const = 1 / 4．
3. 我改变了nfft，到只使用小波或时间向量大小的FFT。这样更容易得到相似振幅的结果。同样，我不把源函数除以nfft，使结果相对独立于所选对象dt．

因此，与你的代码相比，最终的代码只做了很小的更改，变成:

clc;清除所有;vel=2000;介质速度% r=1000;波场观测距离%最终小波出现时间，t=r/vel=0.5秒t= 1;%总时间dt=.0025;%时间步长f0=25;N = round(T/dt);%样本数目t0 = 5/(√(2)*pi*f0);%零时移t = dt*(0:N-1); % time vector tau = t-t0+1*dt; % time vector shifted by t0 src = (1 -2* tau.*tau * f0^2 * pi^2).*exp(-tau.^2 * pi^2 * f0^2); % source wavelet %%FT of signal fmax=1/dt; % max frequency freq= linspace(0,fmax,N); % Freq vector: 0 <= f <= Fmax sw=fft(src); % fourier transform of wavelet, S(w) %% 2-D Greens Function w=2*pi*freq; const = -1i/4; H02 = besselh(0,2,w/vel*r); H02(1)=0; GF = const*H02; %% Calculating the Wavefield WF= ifft(GF.*sw); %% Plot % Wavelet subplot(2,2,1); plot(t,src); title('\bf{Source wavelet}'); xlabel('Time(sec)'); ylabel('Amplitude ') % Frequencies subplot(2,2,2); plot(freq,abs(sw)); title('\bf{Source freq spec}'); xlabel('Frequencies(Hz)'); ylabel('Amplitude ') xlim([0 fmax/2]) % Gneens Fkt subplot(2,2,3); plot(freq,abs(GF)); title('\bf{Green Function}'); xlabel('Frequencies(Hz)'); ylabel('Amplitude ') xlim([0 fmax/2]) % Wafefield subplot(2,2,4); hold on plot(t,real(WF)); title('\bf{Trace}'); xlabel('Time(sec)'); ylabel('Amplitude')

假设信号被采样到奈奎斯特频率，这应该会得到很好的结果: