在MATLAB地震波场的解析解

Question

我一直在试图评估一波的解析解在一个旅行均匀,无限媒体。对于一个给定的源S (t)美元,可以将波场计算距离r美元,对于一个给定的速度v美元

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

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

S(\ω)是美元的频域表示S (t)美元

美元H_0 ^ {(2)} (kr)美元是汉克尔函数的零阶和第二种。

和$ k是波数美元(= \压裂{\ω}{v})美元。

我跟着这链接演示代码在Python中,然而,我不能在MATLAB中实现相同。我也怀疑它是如何工作的。傅里叶变换以来积极以及消极频率(例如$ f{马克斯}/ 2 < f < f{马克斯}/ 2美元)但是代码之间的链接需要频率0 f{马克斯}< f <美元计算Hankal函数。

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

Edit1

作为一个小的修改,我迫使汉克尔函数的第一项为零(因为它是包含一个南实部)。现在工作好但它产生的两座山峰。我认为这是类似于频移效应,作为正面和负面频率FFT后观察,但事实并非如此。在这种情况下振幅是相反的。谁能解释一下吗?

Edit2

如果我认为唯一中的第一个高峰的解决方案是正确的,第二是因为一些影响(我不知道)。问题是当我改变dt在10 ^(例如美元dt \{4}[1、2、3、4、...., 10] $)解决方案的形状和幅度的变化。

任何帮助都是感激。

clc;清除所有;关闭所有;% %参数设置或者= 2000;%的速度中r = 1000;将观察到的波场% %距离最后的小波应该出现在时间,t = r /韦尔= 0.5秒t = 1;%总时间dt =。;%时间步f0 = 25;中央频率% % %创建源小波N =圆(T / dt);%没有样本t0 = 5 /(√(2) *π* 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

所以,我知道很久以前有人问这个问题。但在堆栈的精神交流,我会把我的答案也许未来江南电子竞技平台的用户。

绿党(即功能。、分析解决方案)

首先需要注意的是各种解决方案(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)是一个时间分化:t $ $ f (x) = \压裂{\部分}{\部分t} g (x, t)。$ $这是一个重要的源函数的差异!

绿党函数方程

重申是标准方程(1)在以下形式:$ $ \离开(\压裂{1}{c ^ 2 (x)} \压裂{\部分^ 2}{\部分t ^ 2} - \微分算符^ 2 \右)p (x, t) = f (x, t)。{3}$ $ \标签这个方程转换为频域,美元\部分^ 2 / \部分t ^ 2 \ \ω^ 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 /教学/ AME % 2060633 / / 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 implentation

我做了以下更改w.r.t.代码:

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;清除所有;% %参数设置或者= 2000;%的速度中r = 1000;将观察到的波场% %距离最后的小波应该出现在时间,t = r /韦尔= 0.5秒t = 1;% dt =曾总时间;%时间步f0 = 25;中央频率% % %创建源小波N =圆(T / dt);%没有样本t0 = 5 /(√(2) *π* f0);% 0时移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')

提供信号采样的奈奎斯特频率,这应该给好的结果: