FDTD_CPML_3D.zip_FDTD正演模拟_GPR_GPR正演_fdtd_fdtd 3d

% TM simulation for the 3D FDTD method (CPML) %----------------------------------------------------------------------% tic; close all; clear; clc; % MOVIE = VideoWriter('3D CPML Ricker Movie.avi'); %----------------------------------------------------------------------% c0 = 2.99792458e8; % 光速 eps0 =8.854187818e-12; % 真空中的介电常数 mu0 =4*pi*1e-7; % 真空中的磁导率 %----------------------------------------------------------------------% order = 4; % PML层中的指数系数 thickness = 10; % 10比8效果更佳 xmodel=40; ymodel=40; zmodel=40; xdim = xmodel+2*thickness; % x方向（纵向）总长 ydim = ymodel+2*thickness; % y方向（横向）总长 zdim = zmodel+2*thickness; % y方向（横向）总长 %----------------------------------------------------------------------% freq = 9e8; % 主频 xexcitation = xdim/2; % 源加载的位置 yexcitation = ydim/2; % 源加载的位置 zexcitation = zdim/2; % 源加载的位置 %----------------------------------------------------------------------% % 模拟区域物性参数 eps1=2*ones(xdim,ydim,zdim); % 相对介电常数 sigma1=0.1*ones(xdim,ydim,zdim); % 电导率 %0.05在减少噪声方面比0.005更佳 mu1= ones(xdim,ydim,zdim); % 相对磁导率 %----------------------------------------------------------------------% velocityfield=c0./sqrt(eps1(1:xdim,1:ydim,1:zdim).*mu1(1:xdim,1:ydim,1:zdim)); c_min=min(min(min(velocityfield))); %介质中最小波速 c_max=max(max(max(velocityfield))); %介质中最大波速 %----------------------------------------------------------------------% n_tstep = 500; % 时间总步长 ds= 0.01; % 空间间隔3毫米 % dt= ds/(2*c_max); % 时间间隔 dt=1e-11; %----------------------------------------------------------------------% choose_source=2; %选择激励源 %----------------------------------------------------------------------% if choose_source==1 %-------------------------------雷克子波-------------------------------% omega=2*pi*freq;%omega alfa=0.93*omega;%alfa x=0:1:n_tstep-1; t=x*dt; % 时间 ricker_wavelet=(t).*(t).*exp(-alfa*(t)).*sin(omega*(t)); % 正弦调制的雷克子波 ricker_wavelet=ricker_wavelet/max(abs(ricker_wavelet)); % 归一化 %----------------------------------------------------------------------% % 求激励源的最高频率来确定ds nn = 2^nextpow2(n_tstep); % 总时间步n_tstep所最接近的2的幂次数nn freq_spectrum = abs(fftshift(fft(ricker_wavelet,nn))); % 频谱的振幅值 freq_spectrum = freq_spectrum./max(freq_spectrum); % 归一化 freq_spectrum = freq_spectrum(nn/2+1:end); % 取频谱的后半段 f_max_half = 1/(2*dt); % 频谱中频率最大值的一半 df = 2*f_max_half/nn; % 频谱中频率的采样间隔df f = -f_max_half:df:f_max_half-df; % 横坐标频率 f = f(nn/2+1:end); % 取横坐标频率大于0的部分，即后半段 threshold=0.002; % 阈值，找出最后一个频谱振幅值大于0.002的值作为激励源最高频率 fmax = f(find(freq_spectrum>=threshold, 1, 'last' )); % 找出激励源的最高频率 wavelength_min = c_min/(fmax); % 激励源最高频率所对应的最小波长 ds_limit=wavelength_min/12; % ds应小于ds_limit %----------------------------------------------------------------------% end if choose_source==2 %------------------------------哈里斯脉冲------------------------------% % blackharrispulse % 哈里斯窗的导数 % Note that their formulation has been changed here to T = 1.14/fr, % such that fr represents the approximate dominant frequency of the resulting pulse. % compute the Blackman-Harris window as specified in Chen et al. (1997) x=0:1:n_tstep-1; t=x*dt; % 时间 a = [0.35322222 -0.488 0.145 -0.010222222]; %系数 % T = 1.14/freq; T = 1.55/freq; window = zeros(1,n_tstep); for n=0:3 window = window + a(n+1)*cos(2*n*pi*t./T); end window(t>=T) = 0; % for the pulse, approximate the window's derivative（导数） and normalize（归一化） blackharrispulse = [window(2:end) 0] - window(1:end); % 导数|差分 % blackharrispulse = blackharrispulse./max(abs(blackharrispulse)); %归一化 %----------------------------------------------------------------------% % 求激励源的最高频率来确定ds nn = 2^nextpow2(n_tstep); % 总时间步n_tstep所最接近的2的幂次数nn freq_spectrum = abs(fftshift(fft(blackharrispulse,nn))); % 频谱的振幅值 freq_spectrum = freq_spectrum./max(freq_spectrum); % 归一化 freq_spectrum = freq_spectrum(nn/2+1:end); % 取频谱的后半段 f_max_half = 1/(2*dt); % 频谱中频率最大值的一半 df = 2*f_max_half/nn; % 频谱中频率的采样间隔df f = -f_max_half:df:f_max_half-df; % 横坐标频率 f = f(nn/2+1:end); % 取横坐标频率大于0的部分，即后半段 threshold=0.002; % 阈值，找出最后一个频谱振幅值大于0.002的值作为激励源最高频率 fmax = f(find(freq_spectrum>=threshold, 1, 'last' )); % 找出激励源的最高频率 wavelength_min = c_min/(fmax); % 激励源最高频率所对应的最小波长 ds_limit=wavelength_min/12; % ds应小于ds_limit %----------------------------------------------------------------------% end %----------------------------------------------------------------------% % if ds>ds_limit % msgbox( 'ds应小于ds_limit ','出错!'); % end %----------------------------------------------------------------------% xekappa = ones(1,xdim); % grading of kappa on the electric grid in x-direction xhkappa = ones(1,xdim-1); % grading of kappa on the magnetic grid in x-direction ixekappa = ones(1,xdim); % grading of kappa on the electric grid in x-direction ixhkappa = ones(1,xdim-1); % grading of kappa on the magnetic grid in x-direction xesigma = zeros(1,xdim); % grading of sigma on the electric grid in x-direction xhsigma = zeros(1,xdim-1); % grading of sigma of the magnetic grid in x-direction xealpha = zeros(1,xdim); % alpha就是课本的a_w xhalpha = zeros(1,xdim-1); % grading of alpha on the magnetic grid in x-direction yekappa = ones(1,ydim); % grading of kappa on the electric grid in y-direction yhkappa = ones(1,ydim-1); % grading of kappa on the magnetic grid in y-direction iyekappa = ones(1,ydim); % grading of kappa on the electric grid in y-direction iyhkappa = ones(1,ydim-1); % grading of kappa on the magnetic grid in y-direction yesigma = zeros(1,ydim); % grading of sigma on the electric grid in y-direction yhsigma = zeros(1,ydim-1); % grading of sigma on the magnetic grid in y-direction yealpha = zeros(1,ydim); % grading of alpha on the electric grid in y-direction yhalpha = zeros(1,ydim-1); % grading of alpha on the magnetic grid in y-direction zekappa = ones(1,zdim); % grading of kappa on the electric grid in y-direction zhkappa = ones(1,zdim-1); % grading of kappa on the magnetic grid in y-direction izekappa = ones(1,zdim); % grading of kappa on the electric grid in y-direction izhkappa = ones(1,zdim-1); % grading of kappa on the magnetic grid in y-direction zesigma = zeros(1,zdim); % grading of sigma on the electric grid in y-direction zhsigma = zeros(1,zdim-1); % grading of sigma on the magnetic grid in y-direction zealpha = zeros(1,zdim); % grading of alpha on the electric grid in y-direction zhalpha = zeros(1,zdim-1); % grading of alpha on the magnetic grid in y-direction %-------------------------------------------------------------------------------------------------% sigma_max = (order+1)/(150*sqrt(eps1(thickness+1,thickness+1,thickness+1))*pi*ds); % PML层的最大电阻率(在有些参考书上前面有乘以0.7-1.5的系数) kappa_max = 5; % kappa最大值 ,kappa值的引入是为了改善PML对表面波的吸收特性 alpha_max = 0.008; % 0-0.05， alpha值是为了改善PML对低频分量的吸收特性。 %-------------------------------------------------------------------------------------------------% % 1:10 for i=1:thickness xesigma(1,i)=sigma_max*(((thickness+1)-i)/thickness)^order; % 代入后其实就是(11-i)*ds xekappa(1,i)=1+(kappa_max-1)*(((thickness+1)-i)/thickness)^order; xealpha(1,i)=alpha_max*((i)/(thickness)); %简简单单的等比数列 end % 51:60 for i=xdim-th