CN112051610B - 一种矢量场多模式面波频散计算方法及系统 - Google Patents

一种矢量场多模式面波频散计算方法及系统 Download PDF

Info

Publication number
CN112051610B
CN112051610B CN202011129114.8A CN202011129114A CN112051610B CN 112051610 B CN112051610 B CN 112051610B CN 202011129114 A CN202011129114 A CN 202011129114A CN 112051610 B CN112051610 B CN 112051610B
Authority
CN
China
Prior art keywords
frequency
vector
dispersion
spectrum
surface wave
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202011129114.8A
Other languages
English (en)
Other versions
CN112051610A (zh
Inventor
杨春颖
王赟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Geosciences Beijing
Original Assignee
China University of Geosciences Beijing
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Geosciences Beijing filed Critical China University of Geosciences Beijing
Priority to CN202011129114.8A priority Critical patent/CN112051610B/zh
Publication of CN112051610A publication Critical patent/CN112051610A/zh
Application granted granted Critical
Publication of CN112051610B publication Critical patent/CN112051610B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种矢量场多模式面波频散计算方法及系统。该方法包括：根据检波器采集的数据采用标量MUSIC频谱计算方法确定R分量和Z分量的频散谱，并根据R分量和Z分量的频散谱确定频散曲线的频率区间和速度区间；根据检波器采集的数据采用矢量MUSIC频谱计算方法确定R‑Z矢量场的面波频散谱；在频率区间和速度区间范围内确定面波的频率‑相速度曲线；根据检波器采集的数据采用标量MUSIC频谱计算方法确定Love波的频率‑相速度曲线。采用本发明的方法及系统，能够识别多模式面波，提高了频散计算的精度。

Description

一种矢量场多模式面波频散计算方法及系统
技术领域
本发明涉及面波频散计算技术领域，特别是涉及一种矢量场多模式面波频散计算方法及系统。
背景技术
高密度采集系统提供了高质量地震数据。地球物理工作者能够使用这些高质量的数据，使得重新评估地质结构成为可能。精细、宽频带面波层析将是未来的发展方向。其中，Z分量基阶与其他分量相比，具有较高的分辨率，使得基阶面波反演成为主流的地层层析手段。但是，低速层、非均质性造成面波模式跳跃和模式误判，限制了面波层析的应用。而X分量的面波波场，在能量、频带等方面，与Z分量有一定的差异，可以补充Z分量缺失的面波信息。综合应用多分量、矢量面波波场，将为精细层析提供更准确、更高分辨率的多阶面波，进而提高面波层析的精度。
目前，高分辨率线性拉东变换在多道面波采集系统中得以广泛应用，该方法以线性拉东变换为基础，通过稀疏约束，重新求解面波频散谱。与传统线性拉东变换相比，Z分量面波频散的分辨率可以提高50％。而且增强了高阶面波的识别能力。聚束法是在密集采集阵列基础上发展起来的，该方法使用的面波记录时间较短。假设波场为平面波，根据相邻检波点的传播时间差异，平面波出现相位延迟。以此为基础，通过对相关、不相关信号赋予不同的加权，可以弱化非相关信号，增强相干信号的能量。聚束法输出的功率谱可以描述面波的后方位角和慢度，以此为基础，可以将功率谱转为面波周期-相速度曲线。聚束法包括最小方差无畸变(MVDR)和多重信号分类(MUSIC)等。MUSIC方法具有超高的分辨率。然而，MVDR和标量MUSIC方法主要应用到大尺度地球物理领域，对于面波频散计算的精度不高。
发明内容
本发明的目的是提供一种矢量场多模式面波频散计算方法及系统，能够提高频散计算精度。
为实现上述目的，本发明提供了如下方案：
一种矢量场多模式面波频散计算方法，包括：
获取检波器采集的数据；
根据所述检波器采集的数据采用标量MUSIC频谱计算方法确定R分量和Z分量的频散谱，并根据所述R分量和Z分量的频散谱确定频散曲线的频率区间和速度区间；
根据所述检波器采集的数据采用矢量MUSIC频谱计算方法确定R-Z矢量场的面波频散谱；
在所述频率区间和所述速度区间范围内确定面波的频率-相速度曲线；
根据所述检波器采集的数据采用标量MUSIC频谱计算方法确定Love波的频率-相速度曲线。
可选的，所述根据所述检波器采集的数据采用标量MUSIC频谱计算方法确定R分量和Z分量的频散谱，具体包括：
根据如下公式确定频散谱：
Figure BDA0002734551140000021
其中，
Figure BDA0002734551140000022
A(k)＝[1,exp(-jwτ),...,exp(-j(M-1)wτ)]
Figure BDA0002734551140000023
WL HA(kc)＝1
WL＝μ(R1+LI)-1A(kc)
式中，Pmusic(k)为频散谱，k为波数，k＝2πf/v，f为频率，v为波速，y(t)为检波器的输出，xi(t)为第i个检波器接收的震源子波，ni(t)为第i个检波器的随机噪音，A(k)为导向矩阵，Ai(k)为第i个检波器的导向矩阵，H表示共轭转置，w为角频率，τ为时间，M为检波器总数，R1为标量场的协方差矩阵，Us1为R1分解得到的信号的特征向量，Σs1为R1分解得到的信号的对角阵，Σs1的对角线上的元素为信号的特征值，Un1为R1分解得到的噪音的特征向量，Σn1为R1分解得到的噪音的对角阵，Σn1的对角线上的元素为噪音的特征值，WL为全矢量，μ为常数项，L为加载参数，I为单位矩阵，A(kc)为在面波频散谱的极值点kc处的导向矩阵。
可选的，所述根据所述检波器采集的数据采用矢量MUSIC频谱计算方法确定R-Z矢量场的面波频散谱，具体包括：
根据如下公式确定面波频散谱：
Figure BDA0002734551140000031
其中，
Figure BDA0002734551140000032
Figure BDA0002734551140000033
Z(t)＝F(k)S(t)+N(t)
Figure BDA0002734551140000034
uN＝[coskN,sinkN]T
式中，E表示数学期望，P(k)为面波频散谱，F(k)为矢量阵列的方向向量，R2为矢量场的协方差矩阵，Us2为R2+LI分解得到的信号的特征向量，∑s2为R2+LI分解得到的信号的对角阵，Un2为R2+LI分解得到的噪音的特征向量，∑n2为R2+LI分解得到的噪音的对角阵，Z(t)为检波点接收到的波场矢量，S(t)为震源函数，σ2为检波器白噪音功率，N(t)为背景噪音矢量，f(k)为阵列的方向矢量，N为震源总个数，uN为入射角的正弦和余弦，kN为第N个震源入射到检波器的波数。
可选的，所述在所述频率区间和所述速度区间范围内确定面波的频率-相速度曲线，具体包括：
在所述频率区间和所述速度区间范围内，提取所述R-Z矢量场的面波频散谱最大幅值对应的频散曲线，得到面波的频率-相速度曲线。
本发明还提供一种矢量场多模式面波频散计算系统，包括：
数据获取模块，用于获取检波器采集的数据；
标量MUSIC频谱计算模块，用于根据所述检波器采集的数据采用标量MUSIC频谱计算方法确定R分量和Z分量的频散谱，并根据所述R分量和Z分量的频散谱确定频散曲线的频率区间和速度区间；
矢量MUSIC频谱计算模块，用于根据所述检波器采集的数据采用矢量MUSIC频谱计算方法确定R-Z矢量场的面波频散谱；
面波的频率-相速度曲线确定模块，用于在所述频率区间和所述速度区间范围内确定面波的频率-相速度曲线；
Love波的频率-相速度曲线确定模块，用于根据所述检波器采集的数据采用标量MUSIC频谱计算方法确定Love波的频率-相速度曲线。
可选的，所述标量MUSIC频谱计算模块，具体包括：
标量MUSIC频谱计算单元，用于根据如下公式确定频散谱：
Figure BDA0002734551140000041
其中，
Figure BDA0002734551140000042
A(k)＝[1,exp(-jwτ),...,exp(-j(M-1)wτ)]
Figure BDA0002734551140000043
WL HA(kc)＝1
WL＝μ(R1+LI)-1A(kc)
式中，Pmusic(k)为频散谱，k为波数，k＝2πf/v，f为频率，v为波速，y(t)为检波器的输出，xi(t)为第i个检波器接收的震源子波，ni(t)为第i个检波器的随机噪音，Ai(k)为第i个检波器的导向矩阵，H表示共轭转置，w为角频率，τ为时间，M为检波器总数，R1为标量场的协方差矩阵，Us1为R1分解得到的信号的特征向量，Σs1为R1分解得到的信号的对角阵，Σs1的对角线上的元素为信号的特征值，Un1为R1分解得到的噪音的特征向量，Σn1为R1分解得到的噪音的对角阵，Σn1的对角线上的元素为噪音的特征值，WL为全矢量，μ为常数项，L为加载参数，I为单位矩阵，A(kc)为在面波频散谱的极值点kc处的导向矩阵。
可选的，所述矢量MUSIC频谱计算模块，具体包括：
矢量MUSIC频谱计算单元，用于根据如下公式确定面波频散谱：
Figure BDA0002734551140000044
其中，
Figure BDA0002734551140000045
R2＝E[Z(t)Z(t)H]＝F(k)E[S(t)SH(t)]FH(k)+σ2I
Z(t)＝F(k)S(t)+N(t)
Figure BDA0002734551140000051
uN＝[coskN,sinkN]T
式中，E表示数学期望，P(k)为面波频散谱，F(k)为矢量阵列的方向向量，R2为矢量场的协方差矩阵，Us2为R2+LI分解得到的信号的特征向量，∑s2为R2+LI分解得到的信号的对角阵，Un2为R2+LI分解得到的噪音的特征向量，∑n2为R2+LI分解得到的噪音的对角阵，Z(t)为检波点接收到的波场矢量，S(t)为震源函数，σ2为检波器白噪音功率，N(t)为背景噪音矢量，f(k)为阵列的方向矢量，N为震源总个数，uN为入射角的正弦和余弦，kN为第N个震源入射到检波器的波数。
可选的，所述面波的频率-相速度曲线确定模块，具体包括：
面波的频率-相速度曲线确定单元，用于在所述频率区间和所述速度区间范围内，提取所述R-Z矢量场的面波频散谱最大幅值对应的频散曲线，得到面波的频率-相速度曲线。
与现有技术相比，本发明的有益效果是：
本发明提出了一种矢量场多模式面波频散计算方法及系统，根据检波器采集的数据采用标量MUSIC频谱计算方法确定R分量和Z分量的频散谱，并根据R分量和Z分量的频散谱确定频散曲线的频率区间和速度区间；根据检波器采集的数据采用矢量MUSIC频谱计算方法确定R-Z矢量场的面波频散谱；在频率区间和速度区间范围内确定面波的频率-相速度曲线；根据检波器采集的数据采用标量MUSIC频谱计算方法确定Love波的频率-相速度曲线。本发明适用于规则空间采样的观测系统，同时适用于不规则空间采样的观测系统，能够识别多模式面波，提高了频散计算的精度。此外，本发明使用了对角加载技术，使得频散谱计算过程更稳定。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例中矢量场多模式面波频散计算方法流程图；
图2为本发明实施例中三分面波地震记录示意图；
图3为本发明实施例中均匀道间距测线示意图；
图4为本发明实施例中R、Z单分量面波频散谱图；
图5为本发明实施例中矢量MUSIC法获得R-Z矢量场频谱图；
图6为本发明实施例中T分量Love波频谱图；
图7为本发明实施例中矢量场多模式面波频散计算系统结构图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种矢量场多模式面波频散计算方法及系统，能够提高频散计算精度。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
实施例
图1为本发明实施例中矢量场多模式面波频散计算方法流程图，如图1所示，一种矢量场多模式面波频散计算方法，包括：
步骤101：获取检波器采集的数据。
步骤102：根据检波器采集的数据采用标量MUSIC频谱计算方法确定R分量和Z分量的频散谱，并根据R分量和Z分量的频散谱确定频散曲线的频率区间和速度区间。
步骤102，具体包括：
根据如下公式确定频散谱：
Figure BDA0002734551140000061
其中，
Figure BDA0002734551140000071
A(k)＝[1,exp(-jwτ),...,exp(-j(M-1)wτ)]
Figure BDA0002734551140000072
WL HA(kc)＝1
WL＝μ(R1+LI)-1A(kc)
式中，Pmusic(k)为频散谱，k为波数，k＝2πf/v，f为频率，v为波速，y(t)为检波器的输出，xi(t)为第i个检波器接收的震源子波，ni(t)为第i个检波器的随机噪音，A(k)为导向矩阵，Ai(k)为第i个检波器的导向矩阵，H表示共轭转置，AH(k)为导向矩阵的共轭转置矩阵，w为角频率，τ为时间，M为检波器总数，R1为标量场的协方差矩阵，Us1为R1分解得到的信号的特征向量，Σs1为R1分解得到的信号的对角阵，Σs1的对角线上的元素为信号的特征值，Un1为R1分解得到的噪音的特征向量，Σn1为R1分解得到的噪音的对角阵，Σn1的对角线上的元素为噪音的特征值，WL为全矢量，μ为常数项，L为加载参数，I为单位矩阵，A(kc)为在面波频散谱的极值点kc处的导向矩阵。
步骤103：根据检波器采集的数据采用矢量MUSIC频谱计算方法确定R-Z矢量场的面波频散谱。
步骤103，具体包括：
根据如下公式确定面波频散谱：
Figure BDA0002734551140000073
其中，
Figure BDA0002734551140000074
R2＝E[Z(t)Z(t)H]＝F(k)E[S(t)SH(t)]FH(k)+σ2I
Z(t)＝F(k)S(t)+N(t)
Figure BDA0002734551140000075
uN＝[coskN,sinkN]T
式中，E表示数学期望，P(k)为面波频散谱，F(k)为矢量阵列的方向向量，R2为矢量场的协方差矩阵，Us2为R2+LI分解得到的信号的特征向量，∑s2为R2+LI分解得到的信号的对角阵，Un2为R2+LI分解得到的噪音的特征向量，∑n2为R2+LI分解得到的噪音的对角阵，Z(t)为检波点接收到的波场矢量，S(t)为震源函数，σ2为检波器白噪音功率，N(t)为背景噪音矢量，f(k)为阵列的方向矢量，N为震源总个数，uN为入射角的正弦和余弦，kN为第N个震源入射到检波器的波数。
步骤104：在频率区间和速度区间范围内确定面波的频率-相速度曲线。
步骤104，具体包括：
在频率区间和速度区间范围内，提取R-Z矢量场的面波频散谱最大幅值对应的频散曲线，得到面波的频率-相速度曲线。
步骤105：根据检波器采集的数据采用标量MUSIC频谱计算方法确定Love波的频率-相速度曲线。
本发明首先获得原始面波三分量记录的采集参数；采用标量MUSIC频谱计算方法分别获得R、Z分量的频散谱，了解面波频散基本特征；采用矢量MUSIC方法计算R-Z矢量场的面波频散谱；提取R-Z矢量频散谱最大值对应的频散曲线，保存在文档中；将T分量作为矢量MUSIC的一种特殊情况，计算Love波频散，保存在文档中，完成矢量场多模式频散计算。
本发明提供的矢量场多模式面波频散计算过程，具体内容包括：
步骤1：获得原始面波三分量记录的采集参数。
原始三分面波地震记录，如图2所示，每道的采样点数是500、每条测线的检波点数是8，地震记录的时间采样率4毫秒，道间距是10米。图2(a)为R分量地震记录，图2(b)为Z分量地震记录，图2(c)为T分量地震记录。
步骤2：采用标量MUSIC频谱计算方法分别获得R、Z分量的频散谱。
确定R和Z分量上，可识别频散曲线的频率区间、速度区间、和频散曲线的条数。即，矢量MUSIC谱频散曲线拾取的频率区间、速度区间和可拾取的频散曲线数：
首先，图3为均匀道间距测线示意图，黑色三角表示检波器，波的入射角是θ。如图3所示，对于1个震源8个检波器的均匀道间距、单一测线，接收器的输出为：
Figure BDA0002734551140000081
其中，yi(t)是阵列中第i个检波器的输出，yi(t)可以是R分量接收到的输出，可以是Z分量接收到输出。xi(t)是震源子波，ni(t)是随机噪音，k是平面波的波数，Ai(k)是导向矩阵，表达式为：
A(k)＝[1,exp(-jwτ),...,exp(-j(M-1)wτ)]
式中，w是角频率，τ是时间。当输出yi(t)是Z分量时，时间τ是波数k的函数；当输出yi(t)是R分量时，时间τ是入射θ的正弦函数，波数kc的最优加权向量是：
Wopt＝μR-1A(kc)
其中R是输出信号Y(t)的协方差矩阵，表达式是R＝E[y(t)yH(t)]。根据最优权矢量约束条件wHA(kc)＝1，常数项μ满足:
Figure BDA0002734551140000091
协方差矩阵R可分解为两个子空间:R＝UsΣsUs H+UnΣnUn H，第一项是信号子空间，第二项是噪音子空间。Us是信号的特征向量，Σs是信号的对角阵，对角线上的元素为信号的特征值。Un是信号的特征向量，Σn是噪音的对角阵，对角线上的元素为噪音的特征值。根据信号和噪音子空间的正交性，可获得MUSIC算法功率谱(R、Z分量的频散谱)：
Figure BDA0002734551140000092
然后，将对角加载方法应用到MUSIC方法，采用对角加载，可以衰减噪声波数较小的特征值，以便改善频谱失真问题。进行对角线加载，获得新的权矢量WL：
WL＝μ(R+LI)-1A(kc)
其中，I是单位矩阵，L是加载参数，该参数可使用地震记录的协方差记录来表征。因此，计算过程中，可以使用对角加载公式WL＝μ(R+LI)-1A(kc)替换原公式Wopt＝μR-1A(kc)，增加方法的稳定性。
图4(a)为R分量的频散谱，图4(b)为Z分量的频散谱，面波频散曲线的频率区间是5-30Hz，速度区间是200-450m/s，可识别的频散曲线约为5条。
步骤3：采用矢量MUSIC方法计算R-Z矢量场的面波频散谱。
当面波从震源处向外传播，在r-z平面，波场以方位角θ入射到第m个检波点(xm,ym)。该信号的传输矢量是
Figure BDA0002734551140000101
波数矢量为
Figure BDA0002734551140000102
则
Figure BDA0002734551140000103
式中k是标量波数，假设该检波点的振幅为Am(θ)，则第m检波点记录的波场为
Figure BDA0002734551140000104
即阵列的方向矢量为：
Figure BDA0002734551140000105
符号T表示共轭转置。考虑均匀各向同性完全弹性介质，则各个检波点处具有相同的振幅Am(θ)，方向矢量退化为：
Figure BDA0002734551140000106
当观测系统接收K个震源的入射波场时，则：
Figure BDA0002734551140000107
上式中Y(t)＝[Y1(t),Y2(t),Y3(t),...,YM(t)]是t时刻检波器接收到的波场，波场的方向矩阵F(k)＝[f(k1),f(k2),f(k3),...,f(kM)]，S(t)＝[s1(t),s2(t),s3(t),...,sK(t)]是震源函数，N(t)是背景噪音矢量。平面波在介质中传播，波场传播速度为：
Figure BDA0002734551140000108
上式中，Z是介质波阻抗，X(w)是平面波的频率表达式式，波数k＝2π/λ，λ表示波长。震源到检波点的距离是r，角频率为w，
Figure BDA0002734551140000109
和
Figure BDA00027345511400001010
正交坐标系。在r-z平面，面波z(t)在各分量上的波速与平面波之间的关系为：
z(t)＝[vr,vz]T＝[x(t)cosθ,x(t)sinθ]＝x(t)u(t)
x(t)是时间域的平面波，u(t)＝[cosθ,sinθ]。于是，当n个震源入射到第m个检波器时，该检波器接收到的信号zm(t)为：
Figure BDA00027345511400001011
各个检波点接收到的波场矢量为：
Z(t)＝[z1(t),z2(t),z3(t)...,zM(t)]T
Figure BDA00027345511400001012
其中，um＝[cosθn,sinθn]T，sn(t)符号
Figure BDA00027345511400001013
是克罗内克积。定义矢量阵列的方向向量为：
Figure BDA00027345511400001014
于是，
Z(t)＝F(k)S(t)+N(t)
矢量场的协方差矩阵为：
R＝E[Z(t)Z(t)H]＝F(k)E[S(t)SH(t)]FH(k)+σ2I
式中I是噪音的协方差矩阵。对R进行特征值分解：
Figure BDA0002734551140000111
较大的K个特征值是信号，其余M-K个特征值是噪音。因此，K个特征向量构成信号子空间，M-K个特征向量构成噪音子空间。这两个子空间相互正交。因此，信号方向矢量和噪音子空间正交，矢量阵的功率谱为：
Figure BDA0002734551140000112
当检波器较少，道间距较大时，空间谱的分辨能力降低。针对这种情况，可以采用对角加载方法，对协方差矩阵进行修改，修正后的协方差矩阵为：
RL＝R+LI
式中L是加载量。通过上式，对R进行修正，可以增加矢量场的稳定性，本发明对矢量MUSIC协方差矩阵进行了对角加载。
图5是矢量场频散谱。拾取的区间是5-30Hz，200-450m/s，拾取了5条频散曲线。
步骤4：提取R-Z矢量频散谱最大值对应的频散曲线，保存在文档中。
保存矢量场频散谱振幅最大值对应的频率和相速度，即面波的频率—相速度曲线。
步骤5：将T分量数据看作一个标量场，作为标量MUSIC的输入数据，采用标量MUSIC方法计算Love波频散，保存在文档中，包括：
保存Love波频散谱最大值对应的频率和相速度，即Love波的频率—相速度曲线。
如图6所示，图6为T分量Love波频谱和拾取结果。
步骤6：完成矢量场多模式频散计算。
图7为本发明实施例中矢量场多模式面波频散计算系统结构图。如图7所示，一种矢量场多模式面波频散计算系统，包括：
数据获取模块201，用于获取检波器采集的数据。
标量MUSIC频谱计算模块202，用于根据检波器采集的数据采用标量MUSIC频谱计算方法确定R分量和Z分量的频散谱，并根据R分量和Z分量的频散谱确定频散曲线的频率区间和速度区间。
标量MUSIC频谱计算模块202，具体包括：
标量MUSIC频谱计算单元，用于根据如下公式确定频散谱：
Figure BDA0002734551140000121
其中，
Figure BDA0002734551140000122
A(k)＝[1,exp(-jwτ),...,exp(-j(M-1)wτ)]
Figure BDA0002734551140000123
WL HA(kc)＝1
WL＝μ(R1+LI)-1A(kc)
式中，Pmusic(k)为频散谱，k为波数，k＝2πf/v，f为频率，v为波速，y(t)为检波器的输出，xi(t)为第i个检波器接收的震源子波，ni(t)为第i个检波器的随机噪音，Ai(k)为第i个检波器的导向矩阵，H表示共轭转置，AH(k)为导向矩阵的共轭转置矩阵，w为角频率，τ为时间，M为检波器总数，R1为标量场的协方差矩阵，Us1为R1分解得到的信号的特征向量，Σs1为R1分解得到的信号的对角阵，Σs1的对角线上的元素为信号的特征值，Un1为R1分解得到的噪音的特征向量，Σn1为R1分解得到的噪音的对角阵，Σn1的对角线上的元素为噪音的特征值，WL为全矢量，μ为常数项，L为加载参数，I为单位矩阵，A(kc)为在面波频散谱的极值点kc处的导向矩阵。
矢量MUSIC频谱计算模块203，用于根据检波器采集的数据采用矢量MUSIC频谱计算方法确定R-Z矢量场的面波频散谱。
矢量MUSIC频谱计算模块203，具体包括：
矢量MUSIC频谱计算单元，用于根据如下公式确定面波频散谱：
Figure BDA0002734551140000124
其中，
Figure BDA0002734551140000125
R2＝E[Z(t)Z(t)H]＝F(k)E[S(t)SH(t)]FH(k)+σ2I
Z(t)＝F(k)S(t)+N(t)
Figure BDA0002734551140000131
uN＝[coskN,sinkN]T
式中，E表示数学期望，P(k)为面波频散谱，F(k)为矢量阵列的方向向量，R2为矢量场的协方差矩阵，Us2为R2+LI分解得到的信号的特征向量，∑s2为R2+LI分解得到的信号的对角阵，Un2为R2+LI分解得到的噪音的特征向量，∑n2为R2+LI分解得到的噪音的对角阵，Z(t)为检波点接收到的波场矢量，S(t)为震源函数，σ2为检波器白噪音功率，N(t)为背景噪音矢量，f(k)为阵列的方向矢量，N为震源总个数，uN为入射角的正弦和余弦，kN为第N个震源入射到检波器的波数。
面波的频率-相速度曲线确定模块204，用于在频率区间和速度区间范围内确定面波的频率-相速度曲线。
面波的频率-相速度曲线确定模块204，具体包括：
面波的频率-相速度曲线确定单元，用于在频率区间和速度区间范围内，提取R-Z矢量场的面波频散谱最大幅值对应的频散曲线，得到面波的频率-相速度曲线。
Love波的频率-相速度曲线确定模块205，用于根据检波器采集的数据采用标量MUSIC频谱计算方法确定Love波的频率-相速度曲线。
对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本发明提供了一个完整的矢量场多模式面波频散谱计算过程。矢量MUSIC使用了多分量面波信息，R-Z矢量场获得多模式面波频散曲线。此方法不仅适用于规则观测系统，同时适用于不规则空间采样的观测系统。矢量MUSIC使用了对角加载技术，频散谱计算过程更稳定。本发明将T分量当作矢量MUSIC的一种特殊情况，可以获得T分量的Love波频谱。采用三分量数据，利用矢量MUSIC获得面波的矢量频散谱，这种方法可识别多模式面波，提高了常规Z分量频散计算的精度。本发明中的矢量MUSIC面波频散计算方法，也扩展到各类矢量MUSIC改进算法。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上，本说明书内容不应理解为对本发明的限制。

Claims (8)

1.一种矢量场多模式面波频散计算方法，其特征在于，包括：
获取检波器采集的数据；
根据所述检波器采集的数据采用标量MUSIC频谱计算方法确定R分量和Z分量的频散谱，并根据所述R分量和Z分量的频散谱确定频散曲线的频率区间和速度区间；
根据所述检波器采集的数据采用矢量MUSIC频谱计算方法确定R-Z矢量场的面波频散谱；
在所述频率区间和所述速度区间范围内确定面波的频率-相速度曲线；
根据所述检波器采集的数据采用标量MUSIC频谱计算方法确定Love波的频率-相速度曲线。
2.根据权利要求1所述的矢量场多模式面波频散计算方法，其特征在于，所述根据所述检波器采集的数据采用标量MUSIC频谱计算方法确定R分量和Z分量的频散谱，具体包括：
根据如下公式确定频散谱：
Figure FDA0002734551130000011
其中，
Figure FDA0002734551130000012
A(k)＝[1,exp(-jwτ),...,exp(-j(M-1)wτ)]
Figure FDA0002734551130000013
WL HA(kc)＝1
WL＝μ(R1+LI)-1A(kc)
式中，Pmusic(k)为频散谱，k为波数，k＝2πf/v，f为频率，v为波速，y(t)为检波器的输出，xi(t)为第i个检波器接收的震源子波，ni(t)为第i个检波器的随机噪音，A(k)为导向矩阵，Ai(k)为第i个检波器的导向矩阵，H表示共轭转置，w为角频率，τ为时间，M为检波器总数，R1为标量场的协方差矩阵，Us1为R1分解得到的信号的特征向量，Σs1为R1分解得到的信号的对角阵，Σs1的对角线上的元素为信号的特征值，Un1为R1分解得到的噪音的特征向量，Σn1为R1分解得到的噪音的对角阵，Σn1的对角线上的元素为噪音的特征值，WL为全矢量，μ为常数项，L为加载参数，I为单位矩阵，A(kc)为在面波频散谱的极值点kc处的导向矩阵。
3.根据权利要求2所述的矢量场多模式面波频散计算方法，其特征在于，所述根据所述检波器采集的数据采用矢量MUSIC频谱计算方法确定R-Z矢量场的面波频散谱，具体包括：
根据如下公式确定面波频散谱：
Figure FDA0002734551130000021
其中，
Figure FDA0002734551130000022
R2＝E[Z(t)Z(t)H]＝F(k)E[S(t)SH(t)]FH(k)+σ2I
Z(t)＝F(k)S(t)+N(t)
Figure FDA0002734551130000023
uN＝[coskN,sinkN]T
式中，E表示数学期望，P(k)为面波频散谱，F(k)为矢量阵列的方向向量，R2为矢量场的协方差矩阵，Us2为R2+LI分解得到的信号的特征向量，∑s2为R2+LI分解得到的信号的对角阵，Un2为R2+LI分解得到的噪音的特征向量，∑n2为R2+LI分解得到的噪音的对角阵，Z(t)为检波点接收到的波场矢量，S(t)为震源函数，σ2为检波器白噪音功率，N(t)为背景噪音矢量，f(k)为阵列的方向矢量，N为震源总个数，uN为入射角的正弦和余弦，kN为第N个震源入射到检波器的波数。
4.根据权利要求3所述的矢量场多模式面波频散计算方法，其特征在于，所述在所述频率区间和所述速度区间范围内确定面波的频率-相速度曲线，具体包括：
在所述频率区间和所述速度区间范围内，提取所述R-Z矢量场的面波频散谱最大幅值对应的频散曲线，得到面波的频率-相速度曲线。
5.一种矢量场多模式面波频散计算系统，其特征在于，包括：
数据获取模块，用于获取检波器采集的数据；
标量MUSIC频谱计算模块，用于根据所述检波器采集的数据采用标量MUSIC频谱计算方法确定R分量和Z分量的频散谱，并根据所述R分量和Z分量的频散谱确定频散曲线的频率区间和速度区间；
矢量MUSIC频谱计算模块，用于根据所述检波器采集的数据采用矢量MUSIC频谱计算方法确定R-Z矢量场的面波频散谱；
面波的频率-相速度曲线确定模块，用于在所述频率区间和所述速度区间范围内确定面波的频率-相速度曲线；
Love波的频率-相速度曲线确定模块，用于根据所述检波器采集的数据采用标量MUSIC频谱计算方法确定Love波的频率-相速度曲线。
6.根据权利要求5所述的矢量场多模式面波频散计算系统，其特征在于，所述标量MUSIC频谱计算模块，具体包括：
标量MUSIC频谱计算单元，用于根据如下公式确定频散谱：
Figure FDA0002734551130000031
其中，
Figure FDA0002734551130000032
A(k)＝[1,exp(-jwτ),...,exp(-j(M-1)wτ)]
Figure FDA0002734551130000033
WL HA(kc)＝1
WL＝μ(R1+LI)-1A(kc)
式中，Pmusic(k)为频散谱，k为波数，k＝2πf/v，f为频率，v为波速，y(t)为检波器的输出，xi(t)为第i个检波器接收的震源子波，ni(t)为第i个检波器的随机噪音，Ai(k)为第i个检波器的导向矩阵，H表示共轭转置，w为角频率，τ为时间，M为检波器总数，R1为标量场的协方差矩阵，Us1为R1分解得到的信号的特征向量，Σs1为R1分解得到的信号的对角阵，Σs1的对角线上的元素为信号的特征值，Un1为R1分解得到的噪音的特征向量，Σn1为R1分解得到的噪音的对角阵，Σn1的对角线上的元素为噪音的特征值，WL为全矢量，μ为常数项，L为加载参数，I为单位矩阵，A(kc)为在面波频散谱的极值点kc处的导向矩阵。
7.根据权利要求6所述的矢量场多模式面波频散计算系统，其特征在于，所述矢量MUSIC频谱计算模块，具体包括：
矢量MUSIC频谱计算单元，用于根据如下公式确定面波频散谱：
Figure FDA0002734551130000041
其中，
Figure FDA0002734551130000042
R2＝E[Z(t)Z(t)H]＝F(k)E[S(t)SH(t)]FH(k)+σ2I
Z(t)＝F(k)S(t)+N(t)
Figure FDA0002734551130000043
uN＝[coskN,sinkN]T
式中，E表示数学期望，P(k)为面波频散谱，F(k)为矢量阵列的方向向量，R2为矢量场的协方差矩阵，Us2为R2+LI分解得到的信号的特征向量，∑s2为R2+LI分解得到的信号的对角阵，Un2为R2+LI分解得到的噪音的特征向量，∑n2为R2+LI分解得到的噪音的对角阵，Z(t)为检波点接收到的波场矢量，S(t)为震源函数，σ2为检波器白噪音功率，N(t)为背景噪音矢量，f(k)为阵列的方向矢量，N为震源总个数，uN为入射角的正弦和余弦，kN为第N个震源入射到检波器的波数。
8.根据权利要求7所述的矢量场多模式面波频散计算系统，其特征在于，所述面波的频率-相速度曲线确定模块，具体包括：
面波的频率-相速度曲线确定单元，用于在所述频率区间和所述速度区间范围内，提取所述R-Z矢量场的面波频散谱最大幅值对应的频散曲线，得到面波的频率-相速度曲线。
CN202011129114.8A 2020-10-21 2020-10-21 一种矢量场多模式面波频散计算方法及系统 Active CN112051610B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011129114.8A CN112051610B (zh) 2020-10-21 2020-10-21 一种矢量场多模式面波频散计算方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011129114.8A CN112051610B (zh) 2020-10-21 2020-10-21 一种矢量场多模式面波频散计算方法及系统

Publications (2)

Publication Number Publication Date
CN112051610A CN112051610A (zh) 2020-12-08
CN112051610B true CN112051610B (zh) 2021-06-11

Family

ID=73606530

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011129114.8A Active CN112051610B (zh) 2020-10-21 2020-10-21 一种矢量场多模式面波频散计算方法及系统

Country Status (1)

Country Link
CN (1) CN112051610B (zh)

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6587798B2 (en) * 2000-12-04 2003-07-01 Weatherford/Lamb, Inc. Method and system for determining the speed of sound in a fluid within a conduit
CN104678435A (zh) * 2014-10-27 2015-06-03 李欣欣 一种提取Rayleigh面波频散曲线的方法
CN106646615B (zh) * 2016-12-29 2018-12-25 中国石油天然气集团公司 一种面波频散曲线的数据处理方法及装置
WO2020033465A1 (en) * 2018-08-10 2020-02-13 University Of Houston System Surface wave estimation and removal from seismic data
CN110879410A (zh) * 2019-09-25 2020-03-13 核工业北京地质研究院 一种多分量地震面波勘探方法
CN111103621A (zh) * 2019-12-09 2020-05-05 东华理工大学 一种主动源共成像点叠加多道面波分析方法
CN111736212B (zh) * 2020-07-06 2021-04-20 中国地质大学（北京） 一种假频面波提取方法及系统

Also Published As

Publication number Publication date
CN112051610A (zh) 2020-12-08

Similar Documents

Publication Publication Date Title
Fernandez-Grande et al. A sparse equivalent source method for near-field acoustic holography
US8755249B2 (en) Automatic dispersion extraction of multiple time overlapped acoustic signals
Harley et al. Sparse recovery of the multimodal and dispersive characteristics of Lamb waves
Sun et al. Improving the performance of a vector sensor line array by deconvolution
CN113063490B (zh) 一种基于声压和质点振速双面测量的声场分离方法
Le Courtois et al. Autoregressive model for high-resolution wavenumber estimation in a shallow water environment using a broadband source
CN112285647B (zh) 一种基于稀疏表示与重构的信号方位高分辨估计方法
Ferguson et al. Convolutional neural network for single-sensor acoustic localization of a transiting broadband source in very shallow water
CN109597021B (zh) 一种波达方向估计方法及装置
Salvati et al. Diagonal unloading beamforming in the spherical harmonic domain for acoustic source localization in reverberant environments
Raimondi et al. Tensor decomposition exploiting diversity of propagation velocities: Application to localization of icequake events
CN113608259B (zh) 一种基于iceemdan约束广义s变换的地震薄层检测方法
CN112051610B (zh) 一种矢量场多模式面波频散计算方法及系统
EP0756182A1 (fr) Méthode et dispositif de filtrage d'ondes elliptiques se propageant dans un milieu
CN110261899B (zh) 地震数据z字形干扰波去除方法
CN110703332B (zh) 一种鬼波压制方法
Mahata A subspace algorithm for wideband source localization without narrowband filtering
GB2425355A (en) Residual move out processing of seismic data
CN111736212B (zh) 一种假频面波提取方法及系统
Renfei et al. Method for wavelet denoising of multi-angle prestack seismic data
Soares et al. Environmental inversion using high-resolution matched-field processing
Swanson et al. Small-aperture array processing for passive multi-target angle of arrival estimation
Rozier et al. Inversion of a cylindrical vibrating body in shallow water from aspect-limited data using filtered SVD and the L-curve
Zhai et al. Normal mode energy estimation based on reconstructing the incoherent beamformed outputs from a horizontal array
CN111830565B (zh) 一种基于kl-dsw的tsp多波场分离及噪声压制方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant