提取series中的数值_基于平面波假设的背景噪声数值模拟

背景噪声成像，在有些语境下，也被称为被动源面波成像。该方法实质上是一种通过互相关从弥散的背景噪声场中提取介质的经验格林函数的方法，而经验格林函数中则蕴含着关于介质结构和构成的大量信息。比如，可以通过测量经验格林函数的各个频率的相速度，得到相速度频散曲线。而经典的面波频散分析理论指出，可以通过反演这条频散曲线，得到地下介质横波速度的水平层状结构。为直观地理解这种方法的理论基础，本文基于平面波假设实现某一水平层状模型的背景噪声模拟，并在模拟结果中给出相速度测量值与真实值的对比。

基本原理

我们知道，某一点的时间域波场记录可以理解为震源子波与反射系数的褶积，这在经典的地震勘探理论中，被称为褶积模型。在其他语境中，也被表述为震源子波与格林函数的褶积。时间域的褶积，对应着频率域的乘积。因此，对于记录点A，某一频点的波场记录可以被表示为：

其中，

是子波在该频点的频谱，

是震源位置

到记录点

的格林函数。由平面波假设和面波的几何扩散特征，上式可以具体写为：

"sqrt"代表根号，

是圆波数,可以通过求解层状模型的面波相速度得到；

是依赖于频率，描述介质衰减信息的参数。注意到，在某些频率面波可能有很多相速度，对应于很多阶，此时只需更改e指数部分的波数

，即通过累加的方式得到某一频率的不同阶面波的合成频谱。求得所有频率点的面波频谱信息，对其作傅立叶反变换，就得到了时间域的波形。

算例

背景噪声数值模拟中的理论地质模型被设为两层水平层状模型，模型的参数如表1所示，速度的单位为m/s，密度的单位是

,层厚度的单位是m，在模拟中只考虑基阶波。为尽量接近完全弹性介质，地层衰减因子Q设置为10000（Q越大衰减越小）。需要指出的是，此时得到背景噪声场是瑞雷波记录，分量为Z分量。震源与检波器排列示意图如图1所示，检波器排列呈线性，道间距为150 m,共布设25道。400个震源呈圆环状均匀分布在检波器排列的周围，震源信号被选为雷克子波，子波的延续时间为60 秒，主频随机分布在1~20 Hz之间，最大振幅随机分布在0.001~1之间，子波的起跳时间随机分布在0~60秒之间。对以上参数仿真50次（相当于在实际探测中采集50分钟的记录），得到最后的背景噪声记录，并对其作频散分析。震源主频、最大振幅和起跳时间的分布直方图如图2，图3，图4所示。

表 1 两层水平层状模型参数

横波速度 1200 3000

纵波速度 2400 5200

密度 2000 2300

层厚度 500 半空间

衰减因子 10000

层序号	横波速度	纵波速度	密度	层厚度	衰减因子Q

图 1 震源分布与检波器排列的示意图

图 2 震源主频分布的直方图

图 3 震源最大振幅的分布的直方图

图 4 震源起跳时间的分布直方图

某一检波点处某一仿真过程中形成的背景噪声信号如图5所示，可以看出信号是十分杂乱无章的，而我们的目的就是利用互相关从这种杂乱无章的信号提取出关于介质的信息。以第一道为参考，互相关的时窗长度设为10 s，相邻时窗间重复5 s。最后得到的互相关函数（经验格林函数）炮集如图6所示，图中的互相关信号的振幅已经被归一化。从中可以清晰地看到面波的同相轴，在同相轴上方存在一些模拟中产生的非相干的假象。应用相移法对互相关函数炮集作频散分析，得到的频散能量图如图7所示，图中白色的点线是模型的瑞雷波基阶频散曲线解析解。我们很容易地可以知道，对背景噪声信号经过互相关和叠加处理，最终可以得到来自介质的频散信息。背景噪声成像方法本质上也是一种以互相关和叠加为核心的一种探测介质结构的方法。但是，此处需要指出的是平面波假设只在远场近似条件下才成立，真实世界中，弥散的背景噪声场都是由点源激发的，是柱面波。此时，更符合真实情况的背景噪声模拟方法是离散波数法.

图 5 某一检波点记录的某一时窗的背景噪声信号

图 6 互相关函数炮集

图 7 频散能量图

代码

main.m,用以模拟两层模型的瑞雷波背景噪声场，同时给出互相关炮集。

%% 背景噪声数值模拟，以两层模型为例

clear all;             % 清除matlab的工作空间中的所有变量

% 定义两层模型的参数
VS=[1200 3000];        % 横波速度
VP = [2400 5200];      % 纵波速度
H=[500];               % 地层厚度
den=[2000 2300];       % 地层密度
Q = 10000;             % 地层的衰减因子，此处设置的较大，以接近完全弹性介质

% 设置观测系统参数，定义检波点位置
num_Receiver = 25;
R_loc = zeros(num_Receiver,2); % 检波点的位置矩阵，第一列代表x，第二列代表y
R_loc(:,1) = -1800:150:1800;

% 定义震源位置和参数
SourceNum = 400;
S_r = 5000+1000*rand(1,SourceNum);
S_theta = 2*pi*rand(1,SourceNum);
S_x = S_r.*cos(S_theta);
S_y = S_r.*sin(S_theta);
S_loc = zeros(SourceNum,2);       % 震源点的位置矩阵
S_loc(:,1) = S_x;
S_loc(:,2) = S_y;
fm = [1 20];                      % 子波的主频范围,位于[1 20] Hz
sid = 'ricker';                   % 子波类型，子波是雷克子波
maxamp = [0.001 1];               % 子波的最大振幅范围，位于[0.001 1]

% 设置信号的采样参数
dt = 0.01;                        % 时间域采样间隔
Fs = 1/dt;                        % 信号的采样率
T = 60;                           % 信号总记录时间，秒，即一次仿真过程中的信号长度
nt = round(T/dt);                 % 信号的总采样点数



t=0:dt:(nt-1)*dt;                 % 背景噪声信号的时间轴向量

f = Fs/2*linspace(0,1,nt/2+1);    % 信号的频率轴向量
fmin = 0;                         % 此处代表只模拟[0,25] Hz的面波信息
fmax = 25;
Index = find(f>fmin);
Ind1 = Index(1)-1;
Index = find(f<fmax);
Ind2 = Index(end)+1;

freq = f(Ind1:Ind2);

% 计算两层模型的频散曲线解析解
pv_freq = calcmulti(freq,VS,H,VP,den);
pv_freq = pv_freq(:,1);                % 只考虑基阶面波

ANRec = zeros(nt,num_Receiver);        % 一次仿真过程中的背景噪声记录
T_corr = 10;                           % 设置经验格林函数的时长，为10 s
nt_corr = round(T_corr/dt);
data = zeros(nt_corr,num_Receiver-1);  % 定义以第1道为参考的经验格林函数炮集
data_temp = zeros(nt_corr,num_Receiver-1); % 叠加过程中的中间变量

simu_num = 50;                         % 仿真次数设为50，即相当于观测50*60=3000 s

% 记录仿真过程中所有震源的主频、最大振幅和起跳时间
fm_series = zeros(simu_num*SourceNum,1);
maxamp_series = zeros(simu_num*SourceNum,1);
delay_T_series = zeros(simu_num*SourceNum,1);

tic
count = 0;
for i=1:simu_num
    ANRec(:,:) = 0;
    for j=1:SourceNum
        % 在主频范围内随机生成子波的主频
        fm_temp = fm(1)+(fm(2)-fm(1))*rand;  
        
        % 随机生成子波的起跳时间
        delay_T_temp = T*rand; 
        
        % 随机生成子波的最大振幅
        maxamp_temp = maxamp(1)+(maxamp(2)-maxamp(1))*rand;
        
        % 记录此震源的主频、最大振幅、起跳时间
        count = count+1;
        fm_series(count) = fm_temp;
        maxamp_series(count) = maxamp_temp;
        delay_T_series(count) = delay_T_temp;
        
        % 计算子波
        swave=calcswavelet(fm_temp,delay_T_temp,dt,nt,sid,maxamp_temp);
        swave_spec = fft(swave)/nt;
        swave_spec = swave_spec(1:nt/2+1);
        swave_spec(2:nt/2) = 2*swave_spec(2:nt/2);
        
        S_Spec = swave_spec(Ind1:Ind2);  % 提取需要计算的频带范围内的子波频谱
        
        % 计算此震源产生的背景噪声信号，并叠加
        ANRec_temp = CalcANofSingleS(R_loc,S_loc(j,:),pv_freq,S_Spec,freq,t, Q);
        ANRec = ANRec+ANRec_temp;
    end
    % 对背景噪声数据以第一道为参考作互相关，并令互相关窗口之间重合一半
    [data_temp, Piece_num] = CalcPieceCorr(ANRec,dt,T_corr,0.5, 1);
    data = data+data_temp;
end
t0 = toc;

str = sprintf('主循环计算时间为%f秒.',t0);
disp(str);

% 对经验格林函数炮集按道作归一化
for j=1:num_Receiver-1
    data(:,j) = data(:,j)/max(abs(data(:,j)));
end

% 绘制经验格林函数炮集
[xx,tt] = meshgrid(R_loc(2:end,1),t(1:nt_corr));
figure(1)
pcolor(xx,tt,data);
shading interp;
colormap(jet);
colorbar;
xlabel('Location (m)');
ylabel('Time (s)');
titlestr = sprintf('The cross-correlation gather with %d stacks',simu_num*Piece_num);
title(titlestr);
set(gca,'YDir','reverse');

CalculateDispersionEnergy.m，使用相移法对互相关炮集作频散分析，得到频散能量图。

%% 计算互相关炮集的频散能量图
load data.mat;

% 计算模型的理论频散曲线
VS=[1200 3000];        % 横波速度
VP = [2400 5200];      % 纵波速度
H=[500];               % 地层厚度
den=[2000 2300];       % 地层密度
f = 1:0.25:10;
pv = calcmulti(f,VS,H,VP,den);
pv = pv(:,1);                % 只考虑基阶面波

% 定义频散能量成像的相关参数
dt = 0.01;
offset = 150;
dx = 150;
vmin = 500;
dv = 1;
vmax = 2500;
fmin = 1;
fmax = 10;

% 用相移法形成互相关炮集的频散能量
[E,freq,v] = PhaseShiftOfSW(data,dt,offset,dx,vmin,dv,vmax,fmin,fmax);

figure(1)
pcolor(freq,v,E);
shading interp;
colormap(jet);colorbar;
hold on;
plot(f,pv,'w.','markersize',10);
xlabel('Frequency (Hz)');
ylabel('Phase velocity (m/s)');

---

PS:完整的代码，已经整理到本人的百度企业云盘中，需要的同学可以从下方链接（有效期7天）下载：

ANSimulation​eyun.baidu.com

---