编辑-地球科学江南体育网页版堆栈交换江南电子竞技平台

你不是登录．您的编辑将被放置在队列中，直到完成同行评议．

我们欢迎编辑，使文章更容易理解和更有价值的读者。因为社区成员会审查编辑，所以请尽量让帖子比你发现它的时候更好，例如，通过修改语法或添加额外的资源和超链接。

地震波场的MATLAB解析解

我一直在试着求波在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')

回答

我知道很久以前有人问过这个问题。但本着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)。数组{}\ \端。\tag{2} $$(1)和(2)的区别是一个时间微分:$$ f(x,t) = \frac{\partial}{\partial t}g(x,t)。这是源函数的一个重要区别!##格林函数方程##标准地以以下形式重申方程(1):$$ \left(\frac{1}{c²(x)}\frac{\partial²}{\partial t²}- \nabla²\right)p(x,t) = f(x,t)。\tag{3a}$$将这个方程转换到频域，其中$\partial^2/\partial t^2 \到-\omega^2$，其中$\omega^2/c^2=k^2$，我们得到:$$ \left(-k^2 -\ nabla^2 \right)p(x，\omega) = f(x，\omega)或$$ \left(k^2 + \nabla^2 \right)p(x，\omega) = -f(x，\omega)， \tag{3c}$$我们可以参考，例如，https://www3.nd.edu/~atassi/Teaching/AME%2060633/Notes/greens.pdf来找到Greens函数(其中$f(x，\omega)=\delta(x-x_s)$， $x_s$为源位置)。##格林函数##我只复制结果 For equations of type (1) we have the Green's functions $G$: $$ \begin{array}{lrl} \text{1-D:} & G & = −\frac{i}{2k}e^{−ik\lvert x−x_s \rvert}, \\ \text{2-D:} & G & = −\frac{i}{4}H_0^{(2)}(k \lvert x−x_s \rvert),\\ \text{3-D:} & G & = \frac{1}{4\pi}\frac{e^{−ik\lvert x−x_s \rvert} }{\lvert x−x_s \rvert}. \end{array} $$ Indeed, $H_0^{(2)}$ is the Hankel function, in MATLAB we simply write that with `besselh(0,2, k*r)`. # The MATLAB implentation # I made the following changes w.r.t. your code: 1. The frequency vector can be defined from $f_1\in(0,\frac{1}{\Delta t})$ rather than how you did it, which is $f_2\in(-\frac{1}{2\Delta t},\frac{1}{2\Delta t})$. In this way, the Hankel function only receives positive frequencies. Moreover, it is exactly the order in which the FFT produces its output. I thus also removed the `fftshift` operations. 2. I changed the Green's function to the 2-D case, i.e., I changed `const=-1i*pi` to `const=-1i/4`. 3. I changed the `nfft`, to just using the FFT of the size of the wavelet or time vector. This made it easier to obtain results of similar amplitudes... Also, I don't divide the source function by `nfft`, to make the results relatively independent of the chosen `dt`. Thus, the final code, having made only small changes compared to yours, becomes: clc; clear all; %% Parameter Setting vel=2000; % velocity of medium r=1000; % distance at which wavefield will be observed % the final wavelet should appear at time, t=r/vel=0.5sec T=1; % Total time dt=.0025; % time step f0=25; % central freq %% Create Source Wavelet N = round(T/dt); % No of samples t0 = 5/(sqrt(2)*pi*f0); % Zero time shift 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') Provided that the signal is sampled to the Nyquist frequency, this should give good results: [![Matlab figure results][1]][1] [1]: https://i.stack.imgur.com/ES4um.png

编辑总结

取消

添加评论 |

纠正小的拼写错误
阐明意思而不改变它
添加相关资源或链接
总是尊重作者的意图
不要用编辑回复作者

使用反撇号'或波浪号~创建代码围栏
' ' '
像这样
' ' '
添加语言标识符以突出显示代码
”“python
def函数(foo):
打印(foo)
' ' '
在段落之间放置返回
对于换行符，末尾加2个空格
_italic_或* *的* *
通过将>放在行开始来引用
要建立链接(尽可能使用HTTPS)
< https://example.com >
(例子)(https://example.com)
< a href = " https://example.com " > < / >示例
MathJax方程罪\ \ ^ 2θ美元

格式帮助:
回答帮助:

MathJax help»