CN116400406B - 一种基于阵列的被动源多模式面波频散曲线提取方法 - Google Patents

Abstract

本发明公开了一种基于阵列的被动源多模式面波频散曲线提取方法，包括以下步骤：将密集台站阵列划分为次级台站阵列；计算次级台站阵列的台站对噪声互相关函数；根据次级台站阵列的台站对噪声互相关函数，构建互谱密度矩阵；根据互谱密度矩阵，利用加权互相关波束形成法或修正互相关波束形成法确定频散图；根据频散图提取多模式面波频散曲线。本发明利用台站阵列的背景噪声数据，通过较少的记录时间，提取得到足够精度的多模式面波频散曲线，进而在反演时减小反演的非唯一性，得到高分辨率的速度结构，对大尺度上地球的构造运动和内部过程以及小尺度上断层及其潜在的风险估计提供关键信息。

Description

一种基于阵列的被动源多模式面波频散曲线提取方法

技术领域

本发明涉及地震成像领域，具体涉及一种基于阵列的被动源多模式面波频散曲线提取方法。

背景技术

地震方法是人们认识地球内部结构的主要方法，利用地震方法，地球物理学家可将复杂的地震波形转化为地球模型参数，如地球的速度结构。背景噪声成像方法是21世纪以来地震学研究最为重要的突破性进展，与基于天然地震事件的方法相比，地震干涉方法可从背景噪声的随机场中提取面波频散信息，从而反演出速度结构。具有观测时效高、成像分辨率高、重复性好等优点。随着观测领域中密集台站阵列的广泛布设，促进了台站阵列方法的研究和应用。互相关波束形成(CBF)方法就是一种台站阵列与地震干涉相结合提取面波频散曲线的方法。

但是互相关波束形成(CBF)方法仅简单地将基于天然地震事件的处理过程应用在背景噪声数据中，并没有考虑到背景噪声源随方位分布的情况。其得到的波束图和频散图的分辨率较低，尤其是在高频(或面波波长较短)的情况。因此传统的互相关波束形成方法在提取低模式面波频散曲线时具有不错的效果，但在提取高模式面波频散曲线时，由于频散图分辨率较低无法以一个较高的分辨率识别和提取多模式频散曲线。

发明内容

针对现有技术中的上述不足，本发明提供的一种基于阵列的被动源多模式面波频散曲线提取方法，能利用台站阵列记录的短时间的背景噪声数据提取得到足够精度的多模式面波频散曲线。

为了达到上述发明目的，本发明采用的技术方案为：

一种基于阵列的被动源多模式面波频散曲线提取方法，包括以下步骤：

S1、将密集台站阵列划分为次级台站阵列；

S2、计算步骤S1中次级台站阵列的台站对噪声互相关函数；

S3、根据步骤S2中次级台站阵列的台站对噪声互相关函数，构建互谱密度矩阵；

S4、根据步骤S3中的互谱密度矩阵，利用加权互相关波束形成法或修正互相关波束形成法确定频散图；

S5、根据步骤S4中的频散图提取多模式面波频散曲线。

进一步地，步骤S1包括以下分步骤：

S11、根据期望提取面波的最大波长确定次级台站阵列尺寸，表示为：

L≥2λ

其中：L为阵列尺寸，λ为期望提取面波的最大波长；

S12、根据分步骤S11中的次级台站阵列尺寸将密集台站阵列划分为次级台站阵列，并记录次级台站阵列内的台站位置。

进一步地，步骤S2包括以下分步骤：

S21、获取次级台站阵列的台站对的连续背景噪声记录；

S22、确定窗长和步长；

S23、利用分步骤S22中的窗长和步长将分步骤S21中台站对的连续背景噪声记录划分为不同时间段的背景噪声记录；

S24、将分步骤S23中不同时间段的背景噪声进行互相关并对噪声互相关进行叠加，得到台站对噪声互相关函数，表示为：

cor_ij(ω)＝<d^*(x_i,ω)·d(x_j,ω)>

其中：cor_ij(ω)为频域中台站i和台站j的噪声互相关函数，< >为对长时间的整体平均，d(x_i,ω)为台站i的噪声记录的傅里叶谱，*为取共轭，d(x_j,ω)为台站j的噪声记录的傅里叶谱。

进一步地，在步骤S3中，互谱密度矩阵表示为：

其中：C_ij(ω)为互谱密度矩阵，cor₁₁(ω)为频域中台站1和台站1的噪声互相关函数，cor₁₂(ω)为频域中台站1和台站2的噪声互相关函数，cor_1N(ω)为频域中台站1和台站N的噪声互相关函数，N为次级台站阵列中包含的台站总数，cor₂₁(ω)为频域中台站2和台站1的噪声互相关函数，cor₂₂(ω)为频域中台站2和台站2的噪声互相关函数，cor_N1(ω)为频域中台站N和台站1的噪声互相关函数，cor_NN(ω)为频域中台站N和台站N的噪声互相关函数，为频域中台站1和台站2的噪声互相关函数的共轭，/>为频域中台站1和台站N的噪声互相关函数的共轭。

进一步地，在步骤S4中，利用加权互相关波束形成法确定频散图包括以下步骤：

a1、计算台站间的距离和台站间的方位角；

a2、确定搜索角频率和搜索相速度；

a3、将步骤a1中台站间的距离和台站间的方位角、步骤a2中的搜索频率和搜索相速度以及步骤S3中的互谱密度矩阵代入加权互相关波束形成表达式，搜索得到角频率波束图，表示为：

其中：WCBF(v,θ,ω)为在已知搜索角频率ω下，利用WCBF在搜索相速度v、搜索方向θ上的角频率波束图，N为次级台站阵列中包含的台站总数，i为台站i，j为台站j，C_ij(ω)为互谱密度矩阵，ι为虚数单位，k为波数，θ_ij为台站i指向台站j的方位角，r_ij为台站i和台站j之间的距离；

a4、对分步骤a3中的角频率波束图进行方位叠加得到频散图。

进一步地，在步骤S4中，利用修正互相关波束形成法确定频散图包括以下步骤：

b1、计算台站间的距离；

b2、确定搜索角频率和搜索相速度；

b3、将步骤b1中台站间的距离、步骤b2中的搜索频率和搜索相速度以及步骤S3中的互谱密度矩阵代入修正互相关波束形成表达式，搜索得到频散图。

进一步地，在步骤b3中，第一修正互相关波束形成表达式为：

其中：MCBF(v,ω)指在搜索相速度v和搜索角频率ω下利用MCBF搜索得到的频散图，N为次级台站阵列中包含的台站总数，i为台站i，j为台站j，Re为取实部，C_ij(ω)为互谱密度矩阵，ι为虚数单位，k为波数，r_ij为台站i和台站j之间的距离。

进一步地，在步骤b3中，第二修正互相关波束形成表达式为：

其中：MCBF(v,ω)指在搜索相速度v和搜索角频率ω下利用MCBF搜索得到的频散图，N为次级台站阵列中包含的台站总数，i为台站i，j为台站j，Re为取实部，C_ij(ω)为互谱密度矩阵，J₀为第一类0阶贝塞尔函数，r_ij为台站i和台站j之间的距离，k为波数。

本发明的有益效果为：

(1)本发明利用台站阵列的背景噪声数据，在较少的背景噪声数据记录时间内，通过设计的加权相关波束形成(WCBF)可以得到每个角频率下的波束图，再对波束图进行方位平均可以得到频散图，进而提取得到足够精度的多模式面波频散曲线；

(2)本发明利用台站阵列的背景噪声数据，在较少的背景噪声数据记录时间内，通过设计的修正相关波束形成(MCBF)可以直接得到频散图，并提取得到足够精度的多模式面波频散曲线；

(3)本发明能以高分辨率估计噪声源的分布，在反演时减小反演的非唯一性，得到高分辨率的速度结构，对大尺度上地球的构造运动和内部过程以及小尺度上断层及其潜在的风险估计提供关键信息。

附图说明

图1为一种基于阵列的被动源多模式面波频散曲线提取方法流程图；

图2为台站的噪声记录图，a为台站00的噪声记录图，b为台站01的噪声记录图；

图3为台站00和台站01间的噪声互相关函数在时域上5-20Hz的带通滤波图；

图4为次级台站阵列中每两个台站间的噪声互相关函数按台站对间距排列图；

图5为在搜索角频率为5Hz下加权互相关波束形成的波束图；

图6为每个角频率加权互相关波束形成的角频率波束图进行方位叠加得到的频散图；

图7为修正互相关波束形成的频散图；

图8为多模式面波频散曲线图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，一种基于阵列的被动源多模式面波频散曲线提取方法，包括以下步骤S1-S5：

S1、将密集台站阵列划分为次级台站阵列。

在本发明的一个可选实施例中，本发明取密集台站阵列中一部分位置相近的台站构成次级台站阵列。次级台站阵列中可以包含重复的台站。对于二维分布的大尺寸密集台站阵列，一般以一定的小尺寸进行划分，得到次级台站阵列，以达到对整个区域的充分采样，即移动窗方法。

步骤S1包括以下分步骤：

S11、根据期望提取面波的最大波长确定次级台站阵列尺寸，表示为：

L≥2λ

其中：L为阵列尺寸，即阵列中台站对间的最大距离，λ为期望提取面波的最大波长。

S12、根据分步骤S11中的次级台站阵列尺寸将密集台站阵列划分为次级台站阵列，并记录次级台站阵列内的台站位置。

具体地，本发明对需要记录每个次级台站阵列所包含的每个台站的位置，即每个台站的经度和纬度。

S2、计算步骤S1中次级台站阵列的台站对噪声互相关函数。

在本发明的一个可选实施例中，本发明确定窗长和步长，窗长远大于期望面波最大的周期。本发明确定好窗长，将次级台站阵列中台站对的连续背景噪声记录划分为多个时间段，通过叠加所有时间段内台站对的背景噪声互相关作为两个台站间的噪声互相关函数(NCF)。

步骤S2包括以下分步骤：

S21、获取次级台站阵列的台站对的连续背景噪声记录。

S22、确定窗长和步长。

具体地，窗长根据要探测的最长面波的周期确定，基于采样定理，应远大于2倍的周期。步长取窗长的一半，以充分利用噪声数据，并使下一步骤的噪声互相关函数达到最快速的收敛。

S23、利用分步骤S22中的窗长将分步骤S21中台站对的连续背景噪声记录划分为不同时间段的背景噪声记录。

S24、将分步骤S23中不同时间段的背景噪声进行互相关并对噪声互相关进行叠加，得到台站对噪声互相关函数，表示为：

corij(ω)＝<d*(xi,ω)·d(xj,ω)

其中：cor_ij(ω)为频域中台站i和台站j的噪声互相关函数，可通过反傅里叶变换得到时域中的噪声互相关函数(NCF)，表示为cor_ij(t)，< >为对长时间的整体平均，d(x_i,ω)为台站i的噪声记录的傅里叶谱，*为取共轭，d(x_j,ω)为台站j的噪声记录的傅里叶谱。

具体地，本发明在频域计算每个时间段的两个台站背景噪声数据的互相关，并对所有时间段的噪声互相关进行叠加，得到噪声互相关函数。

本发明为了使噪声互相关函数(NCF)可以在较短的连续背景噪声记录中收敛，在每个台站噪声记录频谱计算的过程中进行非线性操作，包括时域归一化(例如one-bit归一化和滑动平均)和频域归一化(例如滑动平均和仅保留相位)。时域归一化指在时域中进行的归一化操作，以消除异常记录的影响(例如仪器的异常记录和大振幅的地震事件)。频域归一化指在频域中进行的归一化操作，将整个频段的振幅归一化到1单位的幅度，以压制高幅度并增强低幅度的频谱信息。

如图2所示，本发明将台站00和台站01的噪声记录按照窗长40s分割，虚线为分割线，计算每段时间内的互相关，将每个时间段内的互相关进行线性叠加，如图3所示，提供了叠加后的噪声互相关函数在时域上5-20Hz的带通滤波图。如图4所示，本发明提供了次级台站阵列中所有台站对按照台站对间距排列的噪声互相关函数图。

S3、根据步骤S2中次级台站阵列的台站对噪声互相关函数，构建互谱密度矩阵。

在本发明的一个可选实施例中，本发明对划分的每个次级台站阵列，利用步骤S2计算得到的噪声互相关函数(NCF)构建互谱密度矩阵。

互谱密度矩阵表示为：

其中：C_ij(ω)为互谱密度矩阵，cor₁₁(ω)为频域中台站1和台站1的噪声互相关函数，cor₁₂(ω)为频域中台站1和台站2的噪声互相关函数，cor_1N(ω)为频域中台站1和台站N的噪声互相关函数，N为次级台站阵列中包含的台站总数，cor₂₁(ω)为频域中台站2和台站1的噪声互相关函数，cor₂₂(ω)为频域中台站2和台站2的噪声互相关函数，cor_N1(ω)为频域中台站N和台站1的噪声互相关函数，cor_NN(ω)为频域中台站N和台站N的噪声互相关函数，为频域中台站1和台站2的噪声互相关函数的共轭，/>为频域中台站1和台站N的噪声互相关函数的共轭。

S4、根据步骤S3中的互谱密度矩阵，利用加权互相关波束形成法或修正互相关波束形成法确定频散图。

在本发明的一个可选实施例中，本发明计算得到台站间的距离和台站间的方位角，并确定搜索角频率和搜索相速度，根据步骤S3中的互谱密度矩阵，能利用加权互相关波束形成(WCBF)或修正互相关波束形成(MCBF)确定频散图。

利用加权互相关波束形成法确定频散图包括以下步骤：

a1、计算台站间的距离和台站间的方位角。

具体地，本发明计算次级台站阵列中每两个台站间的距离和方位角，对于实际分布的次级台站阵列应按照椭球坐标系计算。次级台站阵列中心为次级阵列中所有台站的中心。本发明以次级阵列中两个台站间的实际距离作为台站间的距离，以其中一个台站指向另一个台站的实际方位作为台站间的方位角。台站间在椭球坐标系中的相对方位角并不是相差180°，因此一个台站对需要计算两次台站间的方位角。

a2、确定搜索角频率和搜索相速度。

a3、将步骤a1中台站间的距离和台站间的方位角、步骤a2中的搜索频率和搜索相速度以及步骤S3中的互谱密度矩阵代入加权互相关波束形成表达式，搜索得到角频率波束图，表示为：

其中：WCBF(v,θ,ω)为在已知搜索角频率ω下，利用WCBF在搜索相速度v、搜索方向θ上的角频率波束图，N为次级台站阵列中包含的台站总数，i为台站i，j为台站j，C_ij(ω)为互谱密度矩阵，ι为虚数单位，k为波数，θ_ij为台站i指向台站j的方位角，r_ij为台站i和台站j之间的距离。

具体地，当搜索角频率ω＝ω₀时，WCBF(v,θ,ω₀)表示在搜索角频率ω₀下加权互相关波束形成的波束图，波束图在搜索相速度v和搜索方向θ上的值表示以相速度v和方位θ入射到次级台站阵列的面波幅值。如图5所示，本发明提供了在搜索角频率为5Hz下加权互相关波束形成的波束图。

a4、对分步骤a3中的角频率波束图进行方位叠加得到频散图。

如图6所示，本发明提供了每个角频率加权互相关波束形成的角频率波束图进行方位叠加得到的频散图。

利用修正互相关波束形成法确定频散图包括以下步骤：

b1、计算台站间的距离。

具体地，本发明计算次级台站阵列中每两个台站间的距离，本发明以次级阵列中两个台站间的实际距离作为台站间的距离。

b2、确定搜索角频率和搜索相速度。

b3、将步骤b1中台站间的距离、步骤b2中的搜索频率和搜索相速度以及步骤S3中的互谱密度矩阵代入修正互相关波束形成表达式，搜索得到频散图。

具体地，本发明能将台站间的距离、搜索频率、搜索相速度以及互谱密度矩阵代入修正互相关波束形成表达式，得到频散图，如图7所示。此处修正互相关波束形成表达式包括第一修正互相关波束形成表达式和第二修正互相关波束形成表达式。

第一修正互相关波束形成表达式为：

其中：MCBF(v,ω)指在搜索相速度v和搜索角频率ω下利用MCBF搜索得到的频散图，N为次级台站阵列中包含的台站总数，i为台站i，j为台站j，Re为取实部，C_ij(ω)为互谱密度矩阵，ι为虚数单位，k为波数，k＝ω/v，ω为搜索角频率，ω＝2πf，f为频率，r_ij为台站i和台站j之间的距离。

第二修正互相关波束形成表达式为：

其中：MCBF(v,ω)指在搜索相速度v和搜索角频率ω下利用MCBF搜索得到的频散图，N为次级台站阵列中包含的台站总数，i为台站i，j为台站j，Re为取实部，C_ij(ω)为互谱密度矩阵，J₀为第一类0阶贝塞尔函数，r_ij为台站i和台站j之间的距离，k为波数，k＝ω/v，ω为搜索角频率，ω＝2πf，f为频率。

S5、根据步骤S4a中的频散图或步骤S4b中的频散图，提取多模式面波频散曲线。

在本发明的一个可选实施例中，频散图中的幅值表示次级台站阵列接收到的由长时间噪声源产生的面波幅值的叠加值。本发明通过提取频散图中识别局部的极值，将连续的极值点的横轴坐标即搜索频率和搜索相速度，作为一个模式面波的频散提取出来，并依次提取多个模式面波。不同模式面波通常表现为频散，即不同搜索角频率以不同的搜索相速度传播(连续变化)。因此在频散图中提取多个连续的频散带(相速度随频率连续变化的极值点)即可提取得到多模式频散曲线，如图8所示。

本发明通过提取多模式面波频散曲线，在反演时减小反演的非唯一性，得到高分辨率的速度结构，对大尺度上地球的构造运动和内部过程以及小尺度上断层及其潜在的风险估计提供关键信息。频散曲线(面波传播相速度随角频率的变化)仅取决于台站阵列下方的介质(速度结构)，而且相同角频率下，不同模式面波的穿透深度是不同的，高模式面波的穿透深度更深。在反演时，不同模式频散是相互独立的，反演的模型必须同时满足多个模式频散曲线，因此使用多个模式可以减小反演的非唯一性。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (7)

1.一种基于阵列的被动源多模式面波频散曲线提取方法，其特征在于，包括以下步骤：

S1、将密集台站阵列划分为次级台站阵列；

S2、计算步骤S1中次级台站阵列的台站对噪声互相关函数；

S3、根据步骤S2中次级台站阵列的台站对噪声互相关函数，构建互谱密度矩阵；

S4、根据步骤S3中的互谱密度矩阵，利用加权互相关波束形成法或修正互相关波束形成法确定频散图；

利用加权互相关波束形成法确定频散图包括以下步骤：

a1、计算台站间的距离和台站间的方位角；

a2、确定搜索角频率和搜索相速度；

a3、将步骤a1中台站间的距离和台站间的方位角、步骤a2中的搜索频率和搜索相速度以及步骤S3中的互谱密度矩阵代入加权互相关波束形成表达式，搜索得到角频率波束图，表示为：

其中：WCBF(v,θ,ω)为在已知搜索角频率ω下，利用WCBF在搜索相速度v、搜索方向θ上的角频率波束图，N为次级台站阵列中包含的台站总数，i为台站i，j为台站j，C_ij(ω)为互谱密度矩阵，ι为虚数单位，k为波数，θ_ij为台站i指向台站j的方位角，r_ij为台站i和台站j之间的距离；

a4、对分步骤a3中的角频率波束图进行方位叠加得到频散图；

S5、根据步骤S4中的频散图提取多模式面波频散曲线。

2.根据权利要求1所述的一种基于阵列的被动源多模式面波频散曲线提取方法，其特征在于，步骤S1包括以下分步骤：

S11、根据期望提取面波的最大波长确定次级台站阵列尺寸，表示为：

L≥2λ

其中：L为阵列尺寸，λ为期望提取面波的最大波长；

S12、根据分步骤S11中的次级台站阵列尺寸将密集台站阵列划分为次级台站阵列，并记录次级台站阵列内的台站位置。

3.根据权利要求1所述的一种基于阵列的被动源多模式面波频散曲线提取方法，其特征在于，步骤S2包括以下分步骤：

S21、获取次级台站阵列的台站对的连续背景噪声记录；

S22、确定窗长和步长；

S23、利用分步骤S22中的窗长和步长将分步骤S21中台站对的连续背景噪声记录划分为不同时间段的背景噪声记录；

S24、将分步骤S23中不同时间段的背景噪声进行互相关并对噪声互相关进行叠加，得到台站对噪声互相关函数，表示为：

cor_ij(ω)＝<d^*(x_i,ω)·d(x_j,ω)>

其中：cor_ij(ω)为频域中台站i和台站j的噪声互相关函数，<>为对长时间的整体平均，d(x_i,ω)为台站i的噪声记录的傅里叶谱，*为取共轭，d(x_j,ω)为台站j的噪声记录的傅里叶谱。

4.根据权利要求1所述的一种基于阵列的被动源多模式面波频散曲线提取方法，其特征在于，在步骤S3中，互谱密度矩阵表示为：

其中：C_ij(ω)为互谱密度矩阵，cor₁₁(ω)为频域中台站1和台站1的噪声互相关函数，cor₁₂(ω)为频域中台站1和台站2的噪声互相关函数，cor_1N(ω)为频域中台站1和台站N的噪声互相关函数，N为次级台站阵列中包含的台站总数，cor₂₁(ω)为频域中台站2和台站1的噪声互相关函数，cor₂₂(ω)为频域中台站2和台站2的噪声互相关函数，cor_N1(ω)为频域中台站N和台站1的噪声互相关函数，cor_NN(ω)为频域中台站N和台站N的噪声互相关函数，为频域中台站1和台站2的噪声互相关函数的共轭，/>为频域中台站1和台站N的噪声互相关函数的共轭。

5.根据权利要求1所述的一种基于阵列的被动源多模式面波频散曲线提取方法，其特征在于，在步骤S4中，利用修正互相关波束形成法确定频散图包括以下步骤：

b1、计算台站间的距离；

b2、确定搜索角频率和搜索相速度；

b3、将步骤b1中台站间的距离、步骤b2中的搜索频率和搜索相速度以及步骤S3中的互谱密度矩阵代入修正互相关波束形成表达式，搜索得到频散图。

6.根据权利要求5所述的一种基于阵列的被动源多模式面波频散曲线提取方法，其特征在于，在步骤b3中，第一修正互相关波束形成表达式为：

其中：MCBF(v,ω)指在搜索相速度v和搜索角频率ω下利用MCBF搜索得到的频散图，N为次级台站阵列中包含的台站总数，i为台站i，j为台站j，Re为取实部，C_ij(ω)为互谱密度矩阵，ι为虚数单位，k为波数，r_ij为台站i和台站j之间的距离。

7.根据权利要求5所述的一种基于阵列的被动源多模式面波频散曲线提取方法，其特征在于，在步骤b3中，第二修正互相关波束形成表达式为：

其中：MCBF(v,ω)指在搜索相速度v和搜索角频率ω下利用MCBF搜索得到的频散图，N为次级台站阵列中包含的台站总数，i为台站i，j为台站j，Re为取实部，C_ij(ω)为互谱密度矩阵，J₀为第一类0阶贝塞尔函数，r_ij为台站i和台站j之间的距离，k为波数。