CN115236738A - 一种基于汉克尔变换的勒夫波频散提取方法 - Google Patents
一种基于汉克尔变换的勒夫波频散提取方法 Download PDFInfo
- Publication number
- CN115236738A CN115236738A CN202210733332.5A CN202210733332A CN115236738A CN 115236738 A CN115236738 A CN 115236738A CN 202210733332 A CN202210733332 A CN 202210733332A CN 115236738 A CN115236738 A CN 115236738A
- Authority
- CN
- China
- Prior art keywords
- wave
- love
- leff
- dispersion
- spectrum
- 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.)
- Pending
Links
- 239000006185 dispersion Substances 0.000 title claims abstract description 134
- 230000009466 transformation Effects 0.000 title claims abstract description 44
- 238000000605 extraction Methods 0.000 title claims abstract description 14
- 238000001228 spectrum Methods 0.000 claims abstract description 113
- 238000000034 method Methods 0.000 claims abstract description 43
- 230000008569 process Effects 0.000 claims abstract description 7
- 238000001514 detection method Methods 0.000 claims abstract description 6
- 238000004613 tight binding model Methods 0.000 claims description 28
- 230000015572 biosynthetic process Effects 0.000 claims description 15
- 238000010586 diagram Methods 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 5
- 239000011159 matrix material Substances 0.000 claims description 4
- 230000003068 static effect Effects 0.000 claims description 4
- 230000010354 integration Effects 0.000 claims description 3
- 238000011160 research Methods 0.000 claims description 3
- 238000013461 design Methods 0.000 claims description 2
- 238000004458 analytical method Methods 0.000 abstract description 14
- 230000010363 phase shift Effects 0.000 description 11
- 238000012545 processing Methods 0.000 description 11
- 238000006073 displacement reaction Methods 0.000 description 10
- 238000005259 measurement Methods 0.000 description 6
- 238000010606 normalization Methods 0.000 description 6
- 238000006243 chemical reaction Methods 0.000 description 4
- 238000003384 imaging method Methods 0.000 description 4
- 238000005070 sampling Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 101150119040 Nsmf gene Proteins 0.000 description 2
- 230000005284 excitation Effects 0.000 description 2
- 239000011435 rock Substances 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 239000002689 soil Substances 0.000 description 2
- 229910000831 Steel Inorganic materials 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 230000036039 immunity Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 239000011148 porous material Substances 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 239000010959 steel Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
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
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- 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
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- 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
- G01V1/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/624—Reservoir parameters
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/63—Seismic attributes, e.g. amplitude, polarity, instant phase
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)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种基于汉克尔变换的勒夫波频散提取方法，具体涉及地震勘探技术领域。本发明通过对野外采集的勒夫波波场记录进行零阶贝塞尔变换，得到勒夫波的频散能量谱，将零阶贝塞尔函数划分为第一类汉克尔函数和第二类汉克尔函数，通过对勒夫波的频散能量谱进行零阶汉克尔变换且变换过程中忽略第二类汉克尔函数，仅在勒夫波波数正方向上对勒夫波波场进行速度扫描，去除勒夫波频散能量谱中的交叉假频，得到变换后的勒夫波频散能量图，提取勒夫波频散能量图中各极值点的相速度值，得到勒夫波的频散曲线。本发明克服了勒夫波波场记录速度扫描时的交叉假频，具有较高的分辨率和多阶分辨能力，为多分量面波探测提供了新的频散分析方法。
Description
技术领域
本发明属于地震勘探技术领域，具体涉及一种基于汉克尔变换的勒夫波频散提取方法。
背景技术
横波速度及其动力学特征作为岩土基本物理性质之一，与岩土岩性、孔隙度、密度和孔隙流体等参数直接相关。由于面波多道分析技术具有非侵入性、高效、低成本、空间采样广和精度高的特点，在获取浅层地表横波速度结构方面得到了广泛应用。相比于瑞雷波，勒夫波的相速度不受纵波速度的影响，勒夫波的反演参数较少使得其反演过程更加稳定、求解的横波速度模型更加可靠。
高分辨率和多模式面波频散分析作为面波多道分析的关键环节，但是勒夫波波场记录进行速度扫描后时会产生交叉假频，交叉假频的存在直接影响了勒夫波频散能量谱的分辨率和多模式分辨能力。因此，亟需提出一种勒夫波频散提取方法消除假频，获取清晰连续的勒夫波频散能量谱，为多分量面波的探测提供新的频散分析方法。
发明内容
本发明为了克服勒夫波波场记录速度扫描时存在交叉假频的问题，解决了勒夫波频散能量谱精度低、连续性差的问题，提出了一种基于汉克尔变换的勒夫波频散提取方法，该方法具有较高的分辨率和多阶分辨能力，为多分量面波的探测提供新的频散分析方法。
为了实现上述目的，本发明采用如下技术方案：
一种基于汉克尔变换的勒夫波频散提取方法，具体包括如下步骤：
步骤1，从研究区野外采集的多道地震记录中提取勒夫波波场记录，对勒夫波波场记录进行零阶贝塞尔变换，得到勒夫波的频散能量谱，如式(1)所示：
式中，r为炮检距，kL为勒夫波的扫描波数，u(ω,r)为炮检距r处的勒夫波波场记录，I为勒夫波的频散能量谱，J0为第一类零阶贝塞尔函数；
步骤2，将零阶贝塞尔函数划分为第一类汉克尔函数和第二类汉克尔函数两个部分，如式(2)所示：
式中，k为波数；
为第一类汉克尔函数，用于表征勒夫波波动方程中的内行柱面波解；
为第二类汉克尔函数，用于表征勒夫波波动方程中的外行柱面波解；J0为第一类零阶贝塞尔函数；
步骤3，将公式(2)代入公式(1)中，对勒夫波的频散能量谱进行零阶汉克尔变换，零阶汉克尔变换过程中忽略第二类汉克尔函数，仅在勒夫波波数的正方向对勒夫波波场进行速度扫描，去除勒夫波频散能量谱中的交叉假频，得到零阶汉克尔变换后的勒夫波频散能量图，在零阶汉克尔变换后的勒夫波频散能量图中各极值点的位置，提取各极值点的相速度值，如式(3)所示：
式中，D(ω,v)为勒夫波的相速度值；u(ω,r)为炮检距r处的勒夫波波场记录；ω为勒夫波的圆频率，ω＝2πf，f为勒夫波的频率；
步骤4，根据勒夫波频散能量图中各极值点的相速度值，提取勒夫波的频散曲线，分析勒夫波的频散特征。
优选地，所述步骤1中，勒夫波波场记录中包括切向格林函数谱和径向格林函数谱两个分量，如式(4)所示：
式中，u(ω,r)为炮检距r处的勒夫波波场记录，gTT(ω,k)为切向格林函数谱，gRR(ω,k)为径向格林函数谱，k为波数，r为炮检距；
由于勒夫波波场记录中既包含切向格林函数谱和径向格林函数谱，也包含零阶贝塞尔函数和一阶贝塞尔函数，所以将勒夫波波场记录经过零阶贝塞尔变换得到的频散能量谱划分为两部分，分别为：
式中，I1为第一部分勒夫波频散能量谱，I2为第二部分勒夫波频散能量谱；
由于波数k和炮检距r相互独立，交换波数k和炮检距r的积分次序，结合同阶贝塞尔函数的正交性，得到：
根据贝塞尔函数的积分性质，得到：
将公式(9)代入公式(8)中，得到：
根据勒夫波波数k→∞时，切向格林函数谱和径向格林函数谱不仅解耦为块对角矩阵，还收敛至对应的静态解，此时切向格林函数谱和径向格林函数谱收敛为：
式中，aL、Δ均为中间参数，
VS为实际地层速度模型中地表的横波速度，VP为实际地层速度模型中地表的纵波速度，μL为实际地层速度模型中地表的剪切模量；
计算实际地层速度模型中勒夫波的最大有效波数值kmax，如式(13)所示：
式中，h为地层厚度；
以实际地层速度模型中勒夫波的最大有效波数值kmax作为分界点，当勒夫波波数k＞kmax时，将公式(11)和公式(12)代入公式(10)中，得到：
式中，A为计算参数，
由于实际地层速度模型为弹性介质，当勒夫波波数k等于勒夫波波数本征值时，切向格林函数谱趋向于无穷大，而当勒夫波波数k等于瑞雷波波数本征值时，径向格林函数谱gRR(ω,k)趋向于无穷大，剔除掉区间[kL,kmax]中所有面波波数的本征值，得到：
并且，
中
和
式中，n为区间[kL,kmax]内瑞雷波的本征模式数量，m为区间[kL,kmax]内勒夫波的本征模式数量；
基于公式(7)和公式(15)可得，当勒夫波扫描波数kL等于勒夫波波数本征值时，第一部分勒夫波频散能量谱I1为无限大，而第二部分勒夫波频散能量谱I2为有限值，勒夫波的频散能量谱I中的极大值与第一部分勒夫波频散能量谱I1中的极大值相对应，因此，利用第一部分勒夫波频散能量谱中的极大值提取勒夫波频散曲线。
优选地，所述步骤3中，勒夫波频散能量谱中交叉假频的波数与检波距成反比，且具有周期性。
本发明所带来的有益技术效果：
本发明提出了一种基于汉克尔变换的勒夫波频散提取方法，利用汉克尔变换提取勒夫波的频散特征，通过改进零阶汉克尔变换将勒夫波波场记录由时间-空间域转换为频率-波数域，并仅在波数正方向上对勒夫波波场进行速度扫描，去除勒夫波频散能量谱中的交叉假频，得到零阶汉克尔变换后的勒夫波频散能量图，获取清晰且连续的勒夫波频散能量谱，解决了勒夫波频散能量谱精度低、连续性差的问题，具有较高的分辨率和多阶分辨能力，为多分量面波的探测提供新的频散分析方法。
附图说明
图1为本发明汉克尔变换的勒夫波频散能量谱对比图。图1中的黑点均为勒夫波的理论相速度，图1中(a)为采用汉克尔函数变换后的勒夫波频散能量谱，图1中(b)为仅采用第一类汉克尔函数变换后的勒夫波频散能量谱。
图2为不同频率下位移切向分量的零阶汉克尔变换结果与切向格林函数谱的对比图。图2中(a)为频率为10Hz时位移切向分量的零阶汉克尔变换结果与切向格林函数谱的对比图，图2中(b)为频率为20Hz时位移切向分量的零阶汉克尔变换结果与切向格林函数谱的对比图，图2中(c)为频率为30Hz时位移切向分量的零阶汉克尔变换结果与切向格林函数谱的对比图，图2中(d)为频率为50Hz时位移切向分量的零阶汉克尔变换结果与切向格林函数谱的对比图。
图3为不同噪声条件下处理的勒夫波频散能量谱。图3中(a)为含有20％随机噪声声的勒夫波波场记录，图3中(b)为含有50％随机噪声声的勒夫波波场记录，图3中(c)为采用传统相移法对含有20％随机噪声声的勒夫波波场记录进行处理得到的勒夫波频散能量谱，图3中(d)为采用传统相移法对含有50％随机噪声声的勒夫波波场记录进行处理得到的勒夫波频散能量谱，图3中(e)为采用本发明方法对含有20％随机噪声声的勒夫波波场记录进行处理得到的勒夫波频散能量谱，图3中(f)为采用本发明方法对含有50％随机噪声声的勒夫波波场记录进行处理得到的勒夫波频散能量谱。
图4为选取测点处的勒夫波波场记录。图4中(a)为测点L1处的勒夫波波场记录，图4中(b)为测点L2处的勒夫波波场记录。
图5为采用本发明方法与传统相移法对测点L1和测点L2处勒夫波波场记录的处理结果图。图5中(a)为采用传统相移法对测点L1处勒夫波波场记录进行处理得到勒夫波的频率能谱图，图5中(b)为采用传统相移法对测点L2处勒夫波波场记录进行处理得到勒夫波的频率能谱图，图5中(c)为采用本发明方法对测点L1处勒夫波波场记录进行处理得到勒夫波的频率能谱图，图5中(d)为采用本发明方法对测点L2处勒夫波波场记录进行处理得到勒夫波的频率能谱图，图5中(e)为测点L1勒夫波波场记录经各道时间域归一化后采用本发明方法生成的频散能量图，图5中(f)为测点L2勒夫波波场记录经各道时间域归一化后采用本发明方法生成的频散能量图。
图6为对选取测点采用不同处理方法时基阶频散曲线和第一高阶频散曲线的提取结果。图6中(a)为针对测点L1采用两种方法对勒夫波基阶和第一高阶的频散曲线提取结果，图6中(b)为针对测点L2采用两种方法对勒夫波基阶和第一高阶的频散曲线提取结果。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明：
本发明提出了一种基于汉克尔变换的勒夫波频散提取方法，具体包括如下步骤：
步骤1，从研究区野外采集的多道地震记录中提取勒夫波波场记录，对勒夫波波场记录进行零阶贝塞尔变换，得到勒夫波的频散能量谱，如式(1)所示：
式中，r为炮检距，kL为勒夫波的扫描波数，u(ω,r)为炮检距r处的勒夫波波场记录，I为勒夫波的频散能量谱，J0为第一类零阶贝塞尔函数；
由于勒夫波波场记录中包括切向格林函数谱和径向格林函数谱两个分量，如式(4)所示：
式中，u(ω,r)为炮检距r处的勒夫波波场记录，gTT(ω,k)为切向格林函数谱，gRR(ω,k)为径向格林函数谱，k为波数，r为炮检距；
由公式(4)可得，勒夫波波场记录中既包含切向格林函数谱和径向格林函数谱，也包含零阶贝塞尔函数和一阶贝塞尔函数，对勒夫波波场记录进行零阶贝塞尔变换后，将得到勒夫波的频散能量谱划分为两个部分：
式中，I1为第一部分勒夫波频散能量谱，I2为第二部分勒夫波频散能量谱；
由于波数k和炮检距r相互独立，交换波数k和炮检距r的积分次序，结合同阶贝塞尔函数的正交性，得到：
根据贝塞尔函数的积分性质，得到：
不失一般性，忽略掉贝塞尔函数中积分性质k＝kL时
的差异，将公式(9)代入公式(8)中，得到：
考虑到勒夫波波数k很大时，切向格林函数谱和径向格林函数谱不仅解耦为块对角矩阵，还收敛至对应的静态解，即当勒夫波波数k→∞时，切向格林函数谱和径向格林函数谱不仅解耦为块对角矩阵，还收敛至对应的静态解，此时切向格林函数谱和径向格林函数谱收敛为：
式中，aL、Δ均为中间参数，
VS为实际地层速度模型中地表的横波速度，VP为实际地层速度模型中地表的纵波速度，μL为实际地层速度模型中地表的剪切模量；
计算实际地层速度模型中勒夫波的最大有效波数值kmax，如式(13)所示：
式中，h为地层厚度；
以实际地层速度模型中勒夫波的最大有效波数值kmax作为分界点，当勒夫波波数k＞kmax时，将公式(11)和公式(12)代入公式(10)中，得到：
式中，A为计算参数，
由于实际地层速度模型为弹性介质，当勒夫波波数k等于勒夫波波数本征值时，切向格林函数谱趋向于无穷大，而当勒夫波波数k等于瑞雷波波数本征值时，径向格林函数谱gRR(ω,k)趋向于无穷大，剔除掉区间[kL,kmax]中所有面波波数的本征值，得到：
并且，
中
和
式中，n为区间[kL,kmax]内瑞雷波的本征模式数量，m为区间[kL,kmax]内勒夫波的本征模式数量；
基于公式(7)和公式(15)可得，当勒夫波扫描波数kL等于勒夫波波数本征值时，第一部分勒夫波频散能量谱I1为无限大，而第二部分勒夫波频散能量谱I2为有限值，即不同勒夫波扫描波数kL条件下，勒夫波的频散能量谱I中的极大值与第一部分勒夫波频散能量谱I1中的极大值相对应。
因此，第一部分勒夫波频散能量谱I1的极大值即对应勒夫波频散曲线，利用第一部分勒夫波频散能量谱中的极大值提取勒夫波频散曲线。
步骤2，将零阶贝塞尔函数划分为第一类汉克尔函数和第二类汉克尔函数两个部分，如式(2)所示：
式中，k为波数；
为第一类汉克尔函数，用于表征勒夫波波动方程中的内行柱面波解；
为第二类汉克尔函数，用于表征勒夫波波动方程中的外行柱面波解；J0为第一类零阶贝塞尔函数。
由公式(2)可得，第一类汉克尔函数
和第二类汉克尔函数
分别描述了二维波动方程的内行柱面波解和外行柱面波解，所以可以将汉克尔变换理解为在勒夫波运行的正方向(+k)和反方向(-k)上分别对勒夫波波场进行速度扫描。
步骤3，由于利用第一类汉克尔函数对勒夫波波场记录进行速度扫描时会产生交叉假频，并且，交叉假频所对应的波数与检波距成反比且具有一定的周期性，为了避免交叉假频的产生，将公式(2)代入公式(1)中对勒夫波的频散能量谱进行零阶汉克尔变换，零阶汉克尔变换过程中忽略第二类汉克尔函数，仅在勒夫波波数正方向上对勒夫波波场进行速度扫描，去除勒夫波频散能量谱中的交叉假频，得到零阶汉克尔变换后的勒夫波频散能量图，在零阶汉克尔变换后的勒夫波频散能量图中各极值点的位置，提取各极值点的相速度值，如式(3)所示：
式中，D(ω,v)为勒夫波的相速度值；u(ω,r)为炮检距r处的勒夫波波场记录；ω为勒夫波的圆频率，ω＝2πf，f为勒夫波的频率。
步骤4，根据勒夫波频散能量图中各极值点的相速度值，提取勒夫波的频散曲线，分析勒夫波的频散特征。
实施例1
本实施例中，采用离散波数法模拟点源引起的地震波场记录，模拟炮集记录的偏移距设置为5m、排列长度设置为200m、道间距设置为5m，模拟得到勒夫波波场记录，采用汉克尔函数对勒夫波波场记录进行汉克尔变化得到的频散能量谱产生交叉假频，如图1(a)所示，而采用本发明方法仅利用第一类汉克尔函数对勒夫波波场记录进行速度扫描增能消除交叉假频，如图1(b)所示，消除交叉假频后得到的频散能量谱更加清晰且连续，更有利于拾取更高频段的频散曲线。
实施例2
将本发明提出的一种基于汉克尔变换的勒夫波频散提取方法应用于双层层状地层模型中，双层层状地层模型的参数如表1所示。
表1双层层状地层模型的参数
对于粘弹性模型，介质的衰减将会导致格林函数谱中的奇点消失，使得面波在无穷大处的值变为有限值。
设置双层层状地层模型中勒夫波的圆频率ω和扫描波数kL，依次将勒夫波的频率设置为10Hz、20Hz、30Hz和50Hz，得到不同频率下位移切向分量的零阶汉克尔变换与切向格林函数谱的对比结果(两者虚部)，如图2所示，结果显示位移切向分量的零阶汉克尔变换与位移切向分量的切向格林函数谱总体差异较小，并且在极大值处两者基本一致，即位移切向分量的零阶汉克尔变换与位移切向分量的切向格林函数谱都等于勒夫波理论相速度，从而验证了利用零阶汉克尔变换可以有效获取勒夫波频散特征。
由于实际应用过程中，采集到的勒夫波波场记录中含有一定的噪声，影响勒夫波频散能量谱的成像质量，所以需要进行抗噪性测试。以表1中的双层层状地层模型为例，分别将20％的随机噪声和50％的随机噪声加入至模拟得到的勒夫波波场记录中，如图3(a)和图3(b)所示，噪声计算采用如下公式：
Data＝Signal+2×(0.5-RND)×mean(max(Signal))×Noise (16)
式中，Data为噪声数据，Signal为模拟的勒夫波波场数据，Noise为噪声水平，mean(max(Signal))为勒夫波波场数据中各道数据最大值的平均值。
图3(c)为采用传统相移法变换生成的勒夫波频散能量图，图3(e)为采用本发明方法变换生成的勒夫波频散能量图，通过对比可得，由于噪声高频干扰，导致在高频段(>40Hz)几乎无法获取有效的频散信息。相比于含有20％随机噪声声的勒夫波波场记录，含有50％随机噪声的勒夫波波场记录所对应的频散能量谱在高频受到的干扰更多，其有效频带范围更窄。通过对比图3(c)、图3(d)和图3(e)、图3(f)可知，采用本发明方法处理得到的勒夫波频散能量谱其抗噪性能优于采用常规相移法处理得到的勒夫波频散能量谱。
应用实验
在雄安新区南部高阳县孝义河河堤土路上应用本发明提出的一种基于汉克尔变换的勒夫波频散提取方法，取得了理想的勒夫波频散提取效果。
实际采集勒夫波波场记录时，采用单端激振法产生勒夫波，为了加强木桩和地面之间的耦合，将多个排钢贯穿桩体并插入地下，激发时多人站在木桩上以增加配重。木桩长轴方向和锤击方向均与检波器方向一致(同垂直于测线方向)，各测点同方向敲击5次，将水平检波器的主频设置为4Hz、检波器间距dx设置为2m、最小偏移距设置为10m、排列长度设置为90m、采样率设置为0.5ms、采样长度设置为1s。
在测线上选取两个测点L1和L2，测点L1和测点L2之间的间距为120m，测量得到测点L1和测点L2处的勒夫波波场记录，如图4(a)和图4(b)所示，通过观察测点L1和测点L2处的勒夫波波场记录可得，勒夫波较发育，信噪比较高。采用传统相移法对测点L1和测点L2处的勒夫波波场记录进行处理得到勒夫波的频率能谱图，如图5(a)和图5(b)所示，采用本发明方法对测点L1和测点L2处勒夫波波场记录进行处理得到勒夫波的频率能谱图，如图5(c)和图5(d)所示。图5中的黑点为根据相移法生成的频散能量图中每个频点的峰值自动拾取的频散曲线，黑色虚线代表线性排列能分辨的最大波长，即直线v/f＝100m。由于勒夫波波场记录的道距较小，根据Nyquist定理，最大波数为Nyquist波数的两倍，为1/dx，也即能分辨的最小波长等于dx。由此可知，拾取的频散曲线并不会受到空间采样假频的干扰。
需要特别说明，为便于描述本实施例简单地按照速度值的大小依次给定频散曲线的阶数，并不考虑面波频散曲线阶数误标的问题。对比图5(a)和图5(c)可得，采用本发明方法对勒夫波第二高阶模式的成像更加清晰连续，进一步对比图5(b)和图5(d)可得，采用本方法对勒夫波第一和第二高阶模式的成像更加清晰连续。分别提取图5(a)和图5(c)中的基阶频散曲线和第一高阶频散曲线，如图6(a)所示，发现两者虽然存在偏差，但是两者的平均差值为1.547m/s，最大差值为4Hz低频端，其值为21.8m/s，这是由于低频段频散能量旁瓣较大，导致误差偏大。而第一高阶两者的最大偏差为4.9m/s，最大相对误差为2.68％，因此，验证了两种方法对勒夫波基阶和第一高阶的成像结果基本一致。再分别提取图5(b)和图5(d)中的基阶频散曲线和第一高阶频散曲线，如图6(b)所示，得到针对测点L2采用两种方法对勒夫波基阶和第一高阶的频散曲线提取结果基本一致。
一般地，面波频散分析前需要对原始勒夫波波场记录进行时间域归一化，由于相移法只关注相位信息，在对原始勒夫波波场记录做Fourier变换后，其振幅谱并没有参与频散分析计算，因此在利用相移法对时间域地震记录进行频散分析前不必进行时间域归一化，而本发明方法同时综合了原始勒夫波波场记录的振幅谱和相位谱，因此时间域归一化可能会对该方法的频散分析结果产生影响。图5(e)和图5(f)分别代表测点L1和测点L2勒夫波波场记录经各道时间域归一化后采用本发明方法生成的频散能量图。通过与图5(c)和图5(d)对比可以看出，由于时间域归一化改变了原始地震信号的振幅谱，导致L1测点第二高阶和L2测点第一高阶只在较窄的频带范围内能量较强，其频散能量整体连续性相对较差。另外，时间域归一化使得基阶频散能量峰值整体有所偏高。
因此，在实际应用中为保证尽量好的频散分析效果，可以利用本发明方法对原始炮集记录和经时间归一化后的炮集记录分别进行频散分析，挑选较好的频散分析结果用于后续的面波反演过程中。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。
Claims (3)
1.一种基于汉克尔变换的勒夫波频散提取方法，其特征在于，具体包括如下步骤：
步骤1，从研究区野外采集的多道地震记录中提取勒夫波波场记录，对勒夫波波场记录进行零阶贝塞尔变换，得到勒夫波的频散能量谱，如式(1)所示：
式中，r为炮检距，kL为勒夫波的扫描波数，u(ω,r)为炮检距r处的勒夫波波场记录，I为勒夫波的频散能量谱，J0为第一类零阶贝塞尔函数；
步骤2，将零阶贝塞尔函数划分为第一类汉克尔函数和第二类汉克尔函数两个部分，如式(2)所示：
式中，k为波数；
为第一类汉克尔函数，用于表征勒夫波波动方程中的内行柱面波解；
为第二类汉克尔函数，用于表征勒夫波波动方程中的外行柱面波解；J0为第一类零阶贝塞尔函数；
步骤3，将公式(2)代入公式(1)中，对勒夫波的频散能量谱进行零阶汉克尔变换，零阶汉克尔变换过程中忽略第二类汉克尔函数，仅在勒夫波波数的正方向对勒夫波波场进行速度扫描，去除勒夫波频散能量谱中的交叉假频，得到零阶汉克尔变换后的勒夫波频散能量图，在零阶汉克尔变换后的勒夫波频散能量图中各极值点的位置，提取各极值点的相速度值，如式(3)所示：
式中，D(ω,v)为勒夫波的相速度值；u(ω,r)为炮检距r处的勒夫波波场记录；ω为勒夫波的圆频率，ω＝2πf，f为勒夫波的频率；
步骤4，根据勒夫波频散能量图中各极值点的相速度值，提取勒夫波的频散曲线，分析勒夫波的频散特征。
2.根据权利要求1所述的一种基于汉克尔变换的勒夫波频散提取方法，其特征在于，所述步骤1中，勒夫波波场记录中包括切向格林函数谱和径向格林函数谱两个分量，如式(4)所示：
式中，u(ω,r)为炮检距r处的勒夫波波场记录，gTT(ω,k)为切向格林函数谱，gRR(ω,k)为径向格林函数谱，k为波数，r为炮检距；
由于勒夫波波场记录中既包含切向格林函数谱和径向格林函数谱，也包含零阶贝塞尔函数和一阶贝塞尔函数，所以将勒夫波波场记录经过零阶贝塞尔变换得到的频散能量谱划分为两部分，分别为：
式中，I1为第一部分勒夫波频散能量谱，I2为第二部分勒夫波频散能量谱；
由于波数k和炮检距r相互独立，交换波数k和炮检距r的积分次序，结合同阶贝塞尔函数的正交性，得到：
根据贝塞尔函数的积分性质，得到：
将公式(9)代入公式(8)中，得到：
根据勒夫波波数k→∞时，切向格林函数谱和径向格林函数谱不仅解耦为块对角矩阵，还收敛至对应的静态解，此时切向格林函数谱和径向格林函数谱收敛为：
式中，aL、Δ均为中间参数，
VS为实际地层速度模型中地表的横波速度，VP为实际地层速度模型中地表的纵波速度，μL为实际地层速度模型中地表的剪切模量；
计算实际地层速度模型中勒夫波的最大有效波数值kmax，如式(13)所示：
式中，h为地层厚度；
以实际地层速度模型中勒夫波的最大有效波数值kmax作为分界点，当勒夫波波数k＞kmax时，将公式(11)和公式(12)代入公式(10)中，得到：
式中，A为计算参数，
由于实际地层速度模型为弹性介质，当勒夫波波数k等于勒夫波波数本征值时，切向格林函数谱趋向于无穷大，而当勒夫波波数k等于瑞雷波波数本征值时，径向格林函数谱gRR(ω,k)趋向于无穷大，剔除掉区间[kL,kmax]中所有面波波数的本征值，得到：
并且，
中
和
式中，n为区间[kL,kmax]内瑞雷波的本征模式数量，m为区间[kL,kmax]内勒夫波的本征模式数量；
基于公式(7)和公式(15)可得，当勒夫波扫描波数kL等于勒夫波波数本征值时，第一部分勒夫波频散能量谱I1为无限大，而第二部分勒夫波频散能量谱I2为有限值，勒夫波的频散能量谱I中的极大值与第一部分勒夫波频散能量谱I1中的极大值相对应，因此，利用第一部分勒夫波频散能量谱中的极大值提取勒夫波频散曲线。
3.根据权利要求1所述的一种基于汉克尔变换的勒夫波频散提取方法，其特征在于，所述步骤3中，勒夫波频散能量谱中交叉假频的波数与检波距成反比，且具有周期性。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210733332.5A CN115236738A (zh) | 2022-06-27 | 2022-06-27 | 一种基于汉克尔变换的勒夫波频散提取方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210733332.5A CN115236738A (zh) | 2022-06-27 | 2022-06-27 | 一种基于汉克尔变换的勒夫波频散提取方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115236738A true CN115236738A (zh) | 2022-10-25 |
Family
ID=83669423
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210733332.5A Pending CN115236738A (zh) | 2022-06-27 | 2022-06-27 | 一种基于汉克尔变换的勒夫波频散提取方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115236738A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115755176A (zh) * | 2022-11-22 | 2023-03-07 | 南方科技大学 | 利用频率汉克尔变换分离波场的面波勘探方法及相关装置 |
-
2022
- 2022-06-27 CN CN202210733332.5A patent/CN115236738A/zh active Pending
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115755176A (zh) * | 2022-11-22 | 2023-03-07 | 南方科技大学 | 利用频率汉克尔变换分离波场的面波勘探方法及相关装置 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| EP1987374B1 (en) | Vh signal integration measure for seismic data | |
| CN101545983A (zh) | 基于小波变换的多属性分频成像方法 | |
| EP2150841A1 (en) | Seismic attributes for reservoir localization | |
| CN106990402A (zh) | 一种基于波浪理论的导航x波段雷达波群检测方法 | |
| CN107153224A (zh) | 检波器动态性能综合测试与评价方法 | |
| CN115236738A (zh) | 一种基于汉克尔变换的勒夫波频散提取方法 | |
| CN102323618B (zh) | 基于分数阶傅里叶变换的相干噪声抑制方法 | |
| Song et al. | Imaging shallow structure with active-source surface wave signal recorded by distributed acoustic sensing arrays | |
| CN104570090B (zh) | 全波形反演噪音滤波算子的提取及使用其噪音滤波的方法 | |
| CN109541689B (zh) | 一种基于反射波能量特征的介质密实度评价方法 | |
| Jiang et al. | Distributed acoustic sensing for shallow structure imaging using mechanical noise: A case study in Guangzhou, China | |
| CN107179548B (zh) | 一种基于真地表的叠前地震成像方法 | |
| CN114415234B (zh) | 基于主动源面波频散和h/v确定浅地表横波速度的方法 | |
| CN110244383B (zh) | 基于近地表数据的地质岩性综合模型创建方法 | |
| CN108646296A (zh) | 基于自适应谱峭度滤波器的沙漠地震信号噪声消减方法 | |
| Smeltzer et al. | Current mapping from the wave spectrum | |
| Taipodia et al. | Resolution of dispersion image obtained from active MASW survey | |
| CN108761534B (zh) | 陆上地震加速度信号应用新方法 | |
| Yang et al. | Imaging dispersion of MASW data with the rammer source | |
| CN108732619B (zh) | 一种海底地球物理数据采集方法 | |
| Hao et al. | Impact of clutter waves on surface-wave analysis under short alignment and a clutter waves removal method | |
| Yang et al. | A frequency-Hankel transform method to extract multimodal Rayleigh wave dispersion spectra from active and passive source surface wave data | |
| deying et al. | Research on seismic data compensation method based on micro-logging Q estimation: A case study in desert area | |
| CN117368967A (zh) | 一种检测用于近地表吸收衰减计算的纵波的方法 | |
| CN116165710A (zh) | 基于面波与纵波的浅层地震分析方法及装置 |
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 |