CN115755176B - 利用频率汉克尔变换分离波场的面波勘探方法及相关装置 - Google Patents

利用频率汉克尔变换分离波场的面波勘探方法及相关装置 Download PDF

Info

Publication number
CN115755176B
CN115755176B CN202211467007.5A CN202211467007A CN115755176B CN 115755176 B CN115755176 B CN 115755176B CN 202211467007 A CN202211467007 A CN 202211467007A CN 115755176 B CN115755176 B CN 115755176B
Authority
CN
China
Prior art keywords
function
wave
hank
frequency
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.)
Active
Application number
CN202211467007.5A
Other languages
English (en)
Other versions
CN115755176A (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.)
Southwest University of Science and Technology
Original Assignee
Southwest University of Science and Technology
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 Southwest University of Science and Technology filed Critical Southwest University of Science and Technology
Priority to CN202211467007.5A priority Critical patent/CN115755176B/zh
Publication of CN115755176A publication Critical patent/CN115755176A/zh
Application granted granted Critical
Publication of CN115755176B publication Critical patent/CN115755176B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明涉及地质勘探技术领域，具体是涉及利用频率汉克尔变换分离波场的面波勘探方法及相关装置。本发明首先采集被勘探区域的地面震动，然后计算勘探区域地层相应的格林函数，之后对格林函数利用频率‑汉克尔变换分离出发散波频谱，最后利用发散波频谱，得到被勘探处所在的地层结构。本发明采用频率汉克尔变换对波场格林函数进行分离，使得分离出来的发散波频谱不包含汇聚波，从而避免了汇聚波对频散谱的噪声干扰，进而使得最后依据频散谱反演出的地层结构更接近于真实的地层结构，提高本发明反演出的地层结构的精度；并实现了主动源(人工产生震源)和被动源(无需人工产生震源，仅采集环境噪声产生的地面振动)的统一，无需改变观测系统和算法。

Description

利用频率汉克尔变换分离波场的面波勘探方法及相关装置
技术领域
本发明涉及地质勘探技术领域，具体是涉及利用频率汉克尔变换分离波场的面波勘探方法及相关装置。
背景技术
面波勘探是工程物探领域应用最广泛的物探方法，比如瑞雷波(Rayleigh wave)这一种面波，其振幅沿纵向垂直于波传播方向呈指数衰减。瑞雷波在层状介质中相速度(波的相位速度)随频率改变而改变，呈现明显的频散特性。水平层状介质中的瑞雷波实际上是纵波P和横波SV在震源区域内各界面处相互干涉叠加而成。瑞雷波在传播过程中能量、频率和速度的变化特征携带了大量地下地层的信息，呈现出的频散特征，也间接反映了层状结构特征。由此研究低频瑞雷波频散可以解决深部地质构造问题；研究高频的瑞雷波可以解决工程勘察、场地和地基处理评价、障碍物和空洞探测等浅层地质问题。
瑞雷波按震源的不同可分为两类：主动源法和被动源法。其中主动源就是人工激发震源产生面波。被动源面波又可成天然源面波，有时又可简称为微动探测和背景噪声成像，顾名思义这一称呼是相对于主动源面波的，即不需要人工主动制造震源，而是将自然界中的潮汐、风、火山活动等自然现象所产生的各种震动，以及来自于人类的各种活动，如车辆行驶、工厂机械运行、人类走动等产生的各种震动作为震源。
在主动源面波勘探中需要人工激发震源，会对环境造成一定影响，且对场地条件有一定的要求。被动源面波勘探的关键问题之一是如何从采集到的地表微动信号中提取面波频散信息，而提取频散信息的关键在于提取地层结构在震源震动下的响应结果所对应的发散波谱，而现有技术提取的发散波谱因包含了汇聚波，汇聚波的干扰会导致最后反演出的地层结构偏离实际的地层结构。
综上所述，现有技术反演出的地层结构精度较低。
因此，现有技术还有待改进和提高。
发明内容
为解决上述技术问题，本发明提供了利用频率汉克尔变换分离波场的面波勘探方法及相关装置，解决了现有技术反演出的地层结构精度较低的问题。
为实现上述目的，本发明采用了以下技术方案：
第一方面，本发明提供一种利用频率汉克尔变换分离波场的面波勘探方法，其中，包括：
计算被勘探处在面波震动下的响应结果所对应的实测格林函数；
对所述实测格林函数应用频率-汉克尔变换分离出所述实测格林函数所涵盖的实测发散波谱；
依据所述实测发散波谱，得到针对所述被勘探处所在的实测地层结构。
在一种实现方式中，所述计算被勘探处在面波震动下的响应结果所对应的实测格林函数，包括：
依据所述面波，得到所述面波中的主动源，所述主动源为人工制造的震源；
依据所述响应结果，得到所述响应结果中的震动数据，所述震动数据为所述被勘探处在所述主动源的震动下产生的震动信息；
将所述震动数据转换到频域，得到频域的所述震动数据；
将所述主动源的震源时间函数转换到频域，得到频域的所述震源时间函数；
依据频域的所述震动数据和频域的所述震源时间函数，计算实测格林函数。
在一种实现方式中，所述计算被勘探处在面波震动下的响应结果所对应的实测格林函数，包括：
依据所述面波，得到所述面波中的被动源，所述被动源为非人工制造的震源；
计算所述被勘探处的任意两个监测设备采集到的两个震动数据之间的互相关函数，得到响应结果；
对所述响应结果应用希尔波特变换，得到变换之后的所述响应结果；
将变换之后的所述响应结果作为实部、所述响应结果作为虚部，得到实测格林函数。
在一种实现方式中，所述实测格林函数为所述被勘探处在面波中的所述主动源或所述被动源震动下的响应结果所对应的函数，所述对所述实测格林函数应用频率-汉克尔变换分离出所述实测格林函数所涵盖的实测发散波谱，包括：
依据所述实测格林函数，得到与所述实测格林函数所对应的地震波场；
计算所述地震波场的周期与波速之积；
计算所述周期与所述波速之积的倒数与距离的乘积结果，所述距离为所述主动源与所述被勘探处之间的距离，或者，所述距离为相邻监测设备之间的距离，所述监测设备用于监测所述被勘探处的震动数据；
计算所述乘积结果的第一类零阶汉克尔函数和第二类零阶汉克尔函数；
依据所述第一类零阶汉克尔函数、所述第二类零阶汉克尔函数、所述实测格林函数，得到频率矢量波数域格林数谱；
提取所述频率矢量波数域格林数谱所涵盖的实测发散波谱。
在一种实现方式中，所述依据所述第一类零阶汉克尔函数、所述第二类零阶汉克尔函数、所述实测格林函数，得到频率矢量波数域格林数谱，包括：
将所述第一类零阶汉克尔函数与所述第二类零阶汉克尔函数之和乘以二分之一，得到贝塞尔函数；
将所述贝塞尔函数、所述实测格林函数与所述距离的乘积关于所述距离做积分，得到频率矢量波数域格林数谱。
在一种实现方式中，所述依据所述实测发散波谱，得到针对所述被勘探处所在的实测地层结构，包括：
依据所述实测发散波谱，绘制图像式的实测频散谱；
对所述实测频散谱应用训练后的神经网络模型，得到针对所述被勘探处所在的实测地层结构。
在一种实现方式中，所述神经网络模型的训练方式，包括：
计算地层结构模型在面波震动下的响应结果所对应的理论格林函数；
计算所述理论格林函数的核函数；
依据所述核函数的频散特性，绘制理论频散谱；
依据所述地层结构模型和所述理论频散谱，训练所述神经网络模型。
第二方面，本发明实施例还提供一种利用频率汉克尔变换分离波场的面波勘探装置，其中，所述装置包括如下组成部分：
格林函数计算模块，用于计算被勘探处在面波震动下的响应结果所对应的实测格林函数；
波谱分离模块，用于对所述实测格林函数应用频率-汉克尔变换分离出所述实测格林函数所涵盖的实测发散波谱；
地层结构计算模块，用于依据所述实测发散波谱，得到针对所述被勘探处所在的实测地层结构。
第三方面，本发明实施例还提供一种终端设备，其中，所述终端设备包括存储器、处理器及存储在所述存储器中并可在所述处理器上运行的利用频率汉克尔变换分离波场的面波勘探程序，所述处理器执行所述利用频率汉克尔变换分离波场的面波勘探程序时，实现上述所述的利用频率汉克尔变换分离波场的面波勘探方法的步骤。
第四方面，本发明实施例还提供一种计算机可读存储介质，所述计算机可读存储介质上存储有利用频率汉克尔变换分离波场的面波勘探程序，所述利用频率汉克尔变换分离波场的面波勘探程序被处理器执行时，实现上述所述的利用频率汉克尔变换分离波场的面波勘探方法的步骤。
有益效果：本发明首先采集被勘探处在面波震动下的响应结果，然后计算响应结果对应的格林函数，之后对格林函数应用频率-汉克尔变换分离出发散波谱，最后根据实测发散波谱，得到针对被勘探处所在的实测地层结构。由于本发明采用频率-汉克尔变换对格林函数进行分离，使得分离出来的发散波谱不包含汇聚波，从而避免了汇聚波对频散谱的噪声干扰，进而使得最后依据频散谱反演出的地层结构更接近于真实的地层结构，即提高了本发明反演出的地层结构的精度。
附图说明
图1为本发明的整体流程图；
图2为本发明实施例中的等间距布设台站示意图；
图3为本发明实施例中的100道合成主动源地震道数据示意图；
图4为实施例中的基于图3中的检波器采集的数据采用F-J方法得到的E所对应的频散谱；
图5为本发明实施例中的基于图3中的检波器采集的数据采用F-H方法得到的发散波谱E1所对应的频散谱；
图6为本发明实施例中的基于图3中的检波器采集的数据得到的汇聚波谱E2；
图7为本发明实施例中的20道合成主动源地震道数据示意图；
图8为实施例中的基于图7中的检波器采集的数据采用F-J方法得到的E所对应的频散谱；
图9为本发明实施例中的基于图7中的检波器采集的数据采用F-H方法得到的发散波谱E1所对应的频散谱；
图10为本发明实施例中的基于图7中的检波器采集的数据得到的汇聚波谱E2；
图11为本发明实施例中的10道合成主动源地震道数据示意图；
图12为实施例中的基于图11中的检波器采集的数据采用F-J方法得到的E所对应的频散谱；
图13为本发明实施例中的基于图11中的检波器采集的数据采用F-H方法得到的发散波谱E1所对应的频散谱；
图14为本发明实施例中的基于图11中的检波器采集的数据得到的汇聚波谱E2；
图15为本发明实施例中的24道实际主动源地震道数据示意图；
图16为本发明实施例中台站和震源分布示意图；
图17为本发明实施例中的基于图15中的检波器采集的数据得到的发散波谱E1所对应的频散谱；
图18为本发明实施例中的基于图15中的检波器采集的数据得到的E所对应的频散谱；
图19为本发明实施例中的基于图15中的检波器采集的数据采用相移法绘制出的频散图；
图20为本发明实施例中的人工锤击震源形成的主动源地震道数据；
图21为本发明实施例中人工锤击的震源与台站之间的分布示意图；
图22为基于图20中的检波器采集的数据绘制的发散波谱E1所对应的频散谱；
图23为基于图20中的检波器采集的数据绘制的发散波谱E1所对应的频散谱；
图24为基于图20中的检波器采集的数据采用相移法绘制的E所对应的频散谱；
图25为本发明实施例中的震源随机分布示意图；
图26为本发明实施例中的合成噪声示意图；
图27为基于图25中的检波器采集的数据绘制的E1的频散谱；
图28为基于图25中的检波器采集的数据绘制的叠加频散谱；
图29为基于图25中的检波器采集的数据绘制的E的频散谱；
图30为基于图25中的检波器采集的数据绘制的E2的频散谱；
图31为本发明的100个震源的分布示意图；
图32为基于图31中的检波器采集的数据绘制的E1的频散谱；
图33为基于图31中的检波器采集的数据绘制的叠加频散谱；
图34为基于图31中的检波器采集的数据绘制的E的频散谱；
图35为基于图31中的检波器采集的数据绘制的E2的频散谱图；
图36为本发明实施例中反演流程图；
图37为本发明实施例提供的终端设备的内部结构原理框图。
具体实施方式
以下结合实施例和说明书附图，对本发明中的技术方案进行清楚、完整地描述。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
经研究发现，面波勘探是工程物探领域应用最广泛的物探方法，比如瑞雷波(Rayleigh wave)这一种面波，其振幅沿纵向垂直于波传播方向呈指数衰减。瑞雷波在层状介质中相速度(波的相位速度)随频率改变而改变，呈现明显的频散特性。水平层状介质中的瑞雷波实际上是纵波P和横波SV在震源区域内各界面处相互干涉叠加而成。瑞雷波在传播过程中能量、频率和速度的变化特征携带了大量地下地层的信息，呈现出的频散特征，也间接反映了层状结构特征。由此研究低频瑞雷波频散可以解决深部地质构造问题；研究高频的瑞雷波可以解决工程勘察、场地和地基处理评价、障碍物和空洞探测等浅层地质问题。瑞雷波按震源的不同可分为两类：主动源法和被动源法。其中主动源就是人工激发震源产生面波。被动源面波又可成天然源面波，有时又可简称为微动探测和背景噪声成像，顾名思义这一称呼是相对于主动源面波的，即不需要人工主动制造震源，而是将自然界中的潮汐、风、火山活动等自然现象所产生的各种震动，以及来自于人类的各种活动，如车辆行驶、工厂机械运行、人类走动等产生的各种震动作为震源。在主动源面波勘探中需要人工激发震源，会对环境造成一定影响，且对场地条件有一定的要求。被动源面波勘探的关键问题之一是如何从采集到的地表微动信号中提取面波频散信息，而提取频散信息的关键在于提取地层结构在震源震动下的响应结果所对应的发散波谱，而现有技术提取的发散波谱因包含了汇聚波，汇聚波的干扰会导致最后反演出的地层结构偏离实际的地层结构。
为解决上述技术问题，本发明提供了利用频率汉克尔变换分离波场的面波勘探方法及相关装置，解决了现有技术反演出的地层结构精度较低的问题。具体实施时，首先计算被勘探处在面波震动下的响应结果所对应的实测格林函数；然后对实测格林函数应用频率-汉克尔变换分离出实测格林函数所涵盖的实测发散波谱；最后依据实测发散波谱，得到针对被勘探处所在的实测地层结构。本发明能够提高反演出的地层结构的精度。
举例说明，比如在地下五十米处设置一个震源，震源的震动带动地表(被勘探处)震动，采集地表的震动数据(比如震动频率和震幅等响应结果)，计算用于表征震动数据和震源二者之间关联度的格林函数(格林函数在地震学中也叫做地震波场，是一种波)，然后，通过频率-汉克尔变换分离出格林函数所涵盖的发散波谱，最后根据发散波谱就可以判断出被勘探处地下的地层结构了。比如，通过发散波谱可以获知，距离被勘探处地下各个深度的波的传播速度，再根据波的传播速度就可以知道各个深度对应的地质材料，各个深度的地质材料就是地层结构。本实施例计算格林函数无需震动数据所包含的波长参与，因此无论震源是主动源还是被动源(主动源和被动源产生的波长不同)在被勘探处产生的震动数据都有相应的格林函数，因此本实施例使得主被动源勘探方法得到了统一。
示例性方法
本实施例的利用频率汉克尔变换分离波场的面波勘探方法可应用于终端设备中，所述终端设备可为具有数据计算功能的终端产品，比如电脑等。在本实施例中，如图1中所示，所述利用频率汉克尔变换分离波场的面波勘探方法具体包括如下步骤：
S100，计算被勘探处在面波震动下的响应结果所对应的实测格林函数。
本实施例的面波包括主动源和被动源，当在被勘探处的各个监测设备都采集到震动数据时，将各个监测设备的震动数据绘制成波形，统计完整波形的数量，根据完整波形的数量判断是主动源还是被动源，完整波形的数量越多，是主动源的概率就越大。
在一个实施例中，主动源为在被勘探处的地表下放置一个震源，震源的震动带动地表震动，通过地表的震动数据(响应结果)和震源时间函数就可以计算出主动源对应的格林函数Gzz(r,w)，该实施例中，步骤S100包括如下的步骤S101至S105：
S101，依据所述面波，得到所述面波中的主动源，所述主动源为人工制造的震源。
S102，依据所述响应结果，得到所述响应结果中的震动数据uz(r,t)，所述震动数据为所述被勘探处在所述主动源的震动下产生的震动信息。
在水平层状各向同性半无限空间模型中，在柱面坐标中，震源(主动源)为原点，地面观测点r(与震源之间的距离为r)观测到的z分量(竖直方向的分量)震动数据uz(r,t)，uz(r,t)是一种波，该波可以用于表示距离震源r处随时间变化的振幅。
S103，将所述震动数据转换到频域，得到频域的所述震动数据。
S104，将所述主动源的震源时间函数转换到频域，得到频域的所述震源时间函数。
S105，依据频域的所述震动数据和频域的所述震源时间函数，计算实测格林函数。
步骤S103是基于如下原理计算出实测格林函数Gzz(r,w)，格林函数也被称为地震波场：
建立震动数据uz(r,t)与震源时间函数fz(t)的关系式。
uz(r,t)＝fz(t)*gzz(r,t) (1)
其中*表示褶积，fz(t)为z分量震源时间函数(震源时间函数的计算为现有技术)，gzz(r,t)为震源与观测点r间z分量的格林函数(时域的格林函数)。将(1)式变换到频率域，可得：
Uz(r,ω)＝Fz(ω)Gzz(r,ω) (2)
对公式(2)进行变形，得到
Gzz(r,ω)＝Uz(r,ω)/Fz(ω) (3)
本实施例中，震源时间函数计算方式包括将震源附近的地震仪记录的数据作为震源时间函数、利用钻井地层数据合成地震子波、利用Ricker子波近似。
在另一个实施例中，面波为被动源，该实施例中，步骤S100包括如下的步骤S106至S109：
S106，依据所述面波，得到所述面波中的被动源，所述被动源为非人工制造的震源。
S107，计算所述被勘探处的任意两个监测设备采集到的两个震动数据之间的互相关函数
Figure BDA0003957834390000101
得到响应结果(/>
Figure BDA0003957834390000102
就是响应结果,r′为两个监测设备之间的距离)。
本实施例中，监测设备为台站(用于接受地表震动信号的地震仪/检波器，在观测系统中称为台站)。多个台站按照三角形、直线型、圆形、不规则方式等排列就构成了台阵，台阵用于计算频散谱。多个台阵组成的观测系统称为测线。台阵勘查中的勘探点，通常为台阵中心点(测点)。如图2所示为沿道路方向等间距布设的台站。
本实施例中的响应结果为两个台站采集的震动数据(震动数据可以是振幅或震动频率)在z方向上的
Figure BDA0003957834390000103
S108，对所述响应结果应用希尔波特变换，得到变换之后的所述响应结果a:
Figure BDA0003957834390000104
Hilbert为希尔波特变换。
Figure BDA0003957834390000105
Figure BDA0003957834390000106
为频率域两个台站间z方向互相关函数。
S109，将变换之后的所述响应结果a作为实部、所述所述响应结果b作为虚部，得到实测格林函数Gzz(r′,w)：
Gzz(r′,w)＝a+bi
式中，i为虚部单位。
S200，对所述实测格林函数应用频率-汉克尔变换分离出所述实测格林函数所涵盖的实测发散波谱。
在一个实施例中，以实测格林函数Gzz(r,w)为例，说明计算实测发散波谱E1(ω,k)的详细过程，该实施例中，步骤S200包括如下的步骤S201至S207：
S201，依据所述实测格林函数Gzz(r,w)，得到与所述实测格林函数所对应的地震波场。
格林函数是在地震学领域中也被称为地震波场。
S202，计算所述地震波场的周期T与波速之积v。
S203，计算所述周期T与所述波速v之积的倒数与距离r的乘积结果x，所述距离r为所述主动源与所述被勘探处之间的距离。
S204，计算所述乘积结果x的第一类零阶汉克尔函数
Figure BDA0003957834390000111
和第二类零阶汉克尔函数/>
Figure BDA0003957834390000112
在一个实施例中，
Figure BDA0003957834390000113
其中λ为波长，ν为波速，ω为频率。
Figure BDA0003957834390000114
代表向外发散的柱面波，/>
Figure BDA0003957834390000115
代表向内汇聚的柱面波。
S205，将所述第一类零阶汉克尔函数与所述第二类零阶汉克尔函数之和乘以二分之一，得到贝塞尔函数J0(x)：
Figure BDA0003957834390000116
S206，将所述贝塞尔函数、所述实测格林函数与所述距离的乘积关于所述距离做积分，得到频率矢量波数域格林数谱E(w,k)：
Figure BDA0003957834390000117
/>
S207，提取所述频率矢量波数域格林数谱所涵盖的实测发散波谱。
利用Hankel函数将E(w,k)分离为发散波谱E1(ω,k)和汇聚波谱E2(ω,k)。
Figure BDA0003957834390000118
Figure BDA0003957834390000119
发散波谱E1(ω,k)的频率、波数/相速度谱能正确的反应面波频散特征，并有效的提高基阶、高阶模式成像质量。而且能够提高去噪能力和成像质量。
本实施例中，公式(6)是基于如下原理形成的：
首先我们对格林函数式进行频率-贝塞尔变换F-J，首先得到频率、矢量波数域Green数谱E(ω,k)：
Figure BDA0003957834390000121
式(9)为S面积分，当在柱面坐标系时(9)可以改写为
Figure BDA0003957834390000122
根据贝塞尔函数积分表达式J0(kr)：
Figure BDA0003957834390000123
将公式(11)代入到公式(10)得到公式(12)
Figure BDA0003957834390000124
将公式(12)结合公式(5)就得到了公式(6)。
在另一个实施例中，在计算Gzz(r′,w)所涵盖的实测发散波谱E′1(w,k)时，只需要将上述公式(6)和(7)中的r替换成r′即可。
S300，依据所述实测发散波谱，得到针对所述被勘探处所在的实测地层结构。
不同的发散波谱对应不同的地层结构，二者之间存在对应关系，因此可以通过获取震源在地层中震动产生的波对应的发散波谱，勘探出地层结构。
在一个实施例中，步骤S300包括如下的步骤S301和S302：
S301，依据所述实测发散波谱，绘制图像式的实测频散谱。
步骤S200得到的实测发散波谱是数学形式，不能直接应用于勘探地层结构，需要绘制成频散谱才能用于勘探地层结构。
S302，对所述实测频散谱应用训练后的神经网络模型，得到针对所述被勘探处所在的实测地层结构。
在一个实施例中，采用图37的流程出反演出地层结构。
训练之后的神经网络模型记录了频散谱与地层结构之间的一一对应关系，因此通过步骤S100和S200得到的发散波谱可以得到与之对应的地层结构。
在一个实施例中，训练神经网络模型包括如下的步骤S3021至S3024：
S3021，计算地层结构模型在面波震动下的响应结果所对应的理论格林函数。
本实施例的地层结构模型是利用某一地区钻孔数据，在实际地层的速度范围、深度范围和厚度范围内随机生成的地层模型，用于训练神经网络的地层模型不小于50000个。在地层结构模型中设置震源，在地层结构模型的表面采集地表面的震动数据(响应结果)，当震源为主动源时就利用公式(3)计算出理论格林函数。当震源为被动源时就利用公式(4)计算出理论格林函数。
S3022，计算所述理论格林函数的核函数
Figure BDA0003957834390000131
S3023，依据所述核函数的频散特性，绘制理论频散谱。
核函数
Figure BDA0003957834390000132
趋于无穷大并具有频散特性，利用格林函数核函数的频散特性，对频率，波数/相速度进行扫描便可绘制地层面波理论频散谱。
S3024，依据所述地层结构模型和所述理论频散谱，训练所述神经网络模型。
通过步骤S3021至S3023生成了大量的训练样本(理论频散谱和地层模型，这二者之间存在对应关系)，将理论频散谱输入至神经网络模型，计算神经网络模型输出的地层结构与地层模型之间的损失值，依据该损失值调整神经网络模型的参数，即实现了训练神经网络模型。
其中步骤S3021至S3023绘制理论频散谱的原理如下：
根据水平层状均匀介质中各向同性震源，格林函数Gzz(r,ω)可表示为
Figure BDA0003957834390000133
其中ω是角频率，k为波数，J0是第一类0阶贝塞尔函数，
Figure BDA0003957834390000134
为格林函数垂直分量的核函数，核函数可以写成分数形式：
Figure BDA0003957834390000135
H(k,ω)是一个光滑没有奇点的函数，而分母D(k,ω)正比久期函数
Figure BDA0003957834390000136
Figure BDA0003957834390000137
(I为单位矩阵，R是反射系数矩阵，D和U分别表示上行和下行，s、l和f表示震源、半空间界面和自由表面。其中久期函数为零时的根为面波频散点，此时波数k＝kn(kn＝ω/cn,cn时第n阶简正振型的相速度)时，核函数/>
Figure BDA0003957834390000141
趋于无穷大并具有频散特性。利用格林函数核函数的频散特性，对频率，波数/相速度进行扫描便可绘制地层面波理论频散谱。
下面以对比实验的方式，证明本发明预测地层结构的准确性：
实验一：绘制在主动源震动作用下产生的频散谱。
a)合成主动源面波地震
采用表1中的数据利用广义反射-投射方法建立模拟的地层模型，然后在模拟的地层模型内部模拟设置垂直地面的单力点源(主动源)，单力点源的震动在地层模型的表面形成地震波场，在地层模型的表面线性排列检波器，最小偏移距1m，数据长度2s，震源时间函数为Ricker子波，主频为15hz。采用不同道间距(检波器之间的距离)的数据，利用式(7)、(8)、(12)得到频率矢量波数域格林数谱E、发散波谱E1(ω,k)和汇聚波谱E2(ω,k)，然后绘制这三者的频散谱。
表1
地层深度(m) 密度(g/cm3) S波速(m/s) P波速度(m/s)
0-10 1.78 180 1500
10-20 1.85 350 1700
20-40 1.80 250 1600
40～ 1.93 600 2000
100道合成如图3所示的主动源地震道数据(用100个检波器采集到的地层表面的震动数据合成地层表面的震动数据或震动波形)。图4是现有技术采用F-J方法得到的E所对应的频散谱。图5是本发明基于F-H方法得到的发散波谱E1所对应的频散谱，图5中的点为理论频散曲线。图6是汇聚波谱E2对应的频散谱。
将图3的100道变成20道(道间距5m)，合成如图7所示的主动源地震道数据。在图7所涵盖的数据基础上，绘制出了图8、图9、图10。其中图8是现有技术采用F-J方法得到的E所对应的频散谱。图9是本发明基于F-H方法得到的发散波谱E1所对应的频散谱。图10是汇聚波谱E2对应的频散谱。
将图3的100道变成10道(道间距10m)，合成如图11所示的主动源地震道数据。在图11所涵盖的数据基础上，绘制出了图12、图13、图14。其中图12是现有技术采用F-J方法得到的E所对应的频散谱。图13是本发明基于F-H方法得到的发散波谱E1所对应的频散谱。图14是汇聚波谱E2对应的频散谱。
从比较上述三组图，可以发现随着台站数量的减少，F-H发散波E1的频散谱可有效的消除汇聚波E2的干扰成像，得到正确的频散谱，并且成像质量比波场不分离成像E提高了很多。
b)实际主动源面波数据
1、数据采集于中国某长江口江滩，震源采用人工锤击，检波器采用等间距线性排列，道间距1m，最小炮检距10m，共24道，合成如图15所示的主动源地震道数据；数据时间长度1s，采样率为1000Hz。现场干扰少，地层均匀，采用rikcer子波作为震源时间函数。台站和震源采用图16中的分布方式，图16中的星号为震源，三角形为台站。绘制出了图17中的发散波谱E1所对应的频散谱、图18中的E所对应的频散谱，采用相移法绘制出图19所示的频散图。
2、数据采集于某市区，震源采用人工锤击，检波器采用等间距线性排列，道间距1m，最小炮检距8m，共24道，长度1s，合成如图20所示的主动源地震道数据，采样频率2000hz。现场干扰大，震源点地层浅层不是很均匀，采用第一道近似为震源时间函数。台站和震源采用图21中的分布方式，图21中的星号为震源，三角形为台站。绘制出了图22中的发散波谱E1所对应的频散谱、图23中的E所对应的频散谱，采用相移法绘制出图24所示的频散图。
从上述来源于两个不同地区的数据可以发现，在干扰程度不同的情况下，我们采用不同的震源时间函数得到了完整了的Green函数，然后通过F-H成像可以有效消除汇聚波干扰提高成像质量，尤其是高阶模式成像，比波场不分离成像E、相移法成像提高很多。
实验二：绘制在主动源震动作用下产生的频散谱。
A.合成远震噪声数据
采用表1四层地层模型合成的如图26所示的60s远震噪声数据，进行频散成像对比分析。如图25所示，1000个随机震源位于半径200～500m圆环内(震源主频15hz)，震源能量和相位均随机分布，图25中星号为随机分布震源，圆圈为台站，100个台站随机分布于半径100m的圆内。基于图25的台站采集的数据绘制出图27中的E1的频散谱，将图27的频散图叠加模型理论频散曲线就形成了图28的频散谱，基于图25的台站采集的数据同样绘制出图29中的E的频散谱以及图30中的E2的频散谱。
B.合成近震数据
实际背景噪声采集中，采集台阵中总是存在少量震源，为了验证F-H对这种情况的有效性和正确性。仍采用表1的四层模型，台阵由24*24共576个台站，台站间距2m，如图31所示，在台阵内分布100个震源(震源主频15hz)，合成的60s噪声数据，研究台阵内少量震源对成像的影响。在实验的过程中，绘制出图32所示的E1的频散谱、图33所示的叠加模型理论频散曲线(在图32的基础上叠加)、图34所示的E的频散谱、图35所示的E2得到的频散谱。
从A和B两组频谱图可以看到，少量震源位于台阵内时，由于合成震源数量仅为100个，所以频散成像质量并没有大量震源位于台站外时好，但F-H发散波E1的频散谱仍比波场不分离成像E的成像质量高；同时也可以看到汇聚波E2不管是主动源还是被动源中的谱都是干扰。
综上，本发明首先采集被勘探处在面波震动下的响应结果，然后计算响应结果对应的格林函数，之后对格林函数应用频率-汉克尔变换分离出发散波谱，最后实测发散波谱，得到针对被勘探处所在的实测地层结构。由于本发明采用频率-汉克尔变换对格林函数进行分离，使得分离出来的发散波谱不包含汇聚波，从而避免了汇聚波对频散谱的噪声干扰，进而使得最后依据频散谱反演出的地层结构更接近于真实的地层结构，即提高了本发明反演出的地层结构的精度。
示例性装置
本实施例还提供一种利用频率汉克尔变换分离波场的面波勘探装置，所述装置包括如下组成部分：
格林函数计算模块，用于计算被勘探处在面波震动下的响应结果所对应的实测格林函数；
波谱分离模块，用于对所述实测格林函数应用频率-汉克尔变换分离出所述实测格林函数所涵盖的实测发散波谱；
地层结构计算模块，用于依据所述实测发散波谱，得到针对所述被勘探处所在的实测地层结构。
基于上述实施例，本发明还提供了一种终端设备，其原理框图可以如图37所示。该终端设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏、温度传感器。其中，该终端设备的处理器用于提供计算和控制能力。该终端设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该终端设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种利用频率汉克尔变换分离波场的面波勘探方法。该终端设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该终端设备的温度传感器是预先在终端设备内部设置，用于检测内部设备的运行温度。
本领域技术人员可以理解，图37中示出的原理框图，仅仅是与本发明方案相关的部分结构的框图，并不构成对本发明方案所应用于其上的终端设备的限定，具体的终端设备以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
在一个实施例中，提供了一种终端设备，终端设备包括存储器、处理器及存储在存储器中并可在处理器上运行的利用频率汉克尔变换分离波场的面波勘探程序，处理器执行利用频率汉克尔变换分离波场的面波勘探程序时，实现如下操作指令：
计算被勘探处在面波震动下的响应结果所对应的实测格林函数；
对所述实测格林函数应用频率-汉克尔变换分离出所述实测格林函数所涵盖的实测发散波谱；
依据所述实测发散波谱，得到针对所述被勘探处所在的实测地层结构。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本发明所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(ROM)、可编程ROM(PROM)、电可编程ROM(EPROM)、电可擦除可编程ROM(EEPROM)或闪存。易失性存储器可包括随机存取存储器(RAM)或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM(SRAM)、动态RAM(DRAM)、同步DRAM(SDRAM)、双数据率SDRAM(DDRSDRAM)、增强型SDRAM(ESDRAM)、同步链路(Synchlink)DRAM(SLDRAM)、存储器总线(Rambus)直接RAM(RDRAM)、直接存储器总线动态RAM(DRDRAM)、以及存储器总线动态RAM(RDRAM)等。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (9)

1.一种利用频率汉克尔变换分离波场的面波勘探方法，其特征在于，包括：
计算被勘探处在面波震动下的响应结果所对应的实测格林函数；
对所述实测格林函数应用频率-汉克尔变换分离出所述实测格林函数所涵盖的实测发散波谱；
依据所述实测发散波谱，得到针对所述被勘探处所在的实测地层结构；
所述实测格林函数为所述被勘探处在面波中的主动源或被动源震动下的响应结果所对应的函数，所述对所述实测格林函数应用频率-汉克尔变换分离出所述实测格林函数所涵盖的实测发散波谱，包括：
依据所述实测格林函数，得到与所述实测格林函数所对应的地震波场；
计算所述地震波场的周期与波速之积；
计算所述周期与所述波速之积的倒数与距离的乘积结果，所述距离为所述主动源与所述被勘探处之间的距离，或者，所述距离为相邻监测设备之间的距离，所述监测设备用于监测所述被勘探处的震动数据；
计算所述乘积结果的第一类零阶汉克尔函数和第二类零阶汉克尔函数；
依据所述第一类零阶汉克尔函数、所述第二类零阶汉克尔函数、所述实测格林函数，得到频率矢量波数域格林数谱；
提取所述频率矢量波数域格林数谱所涵盖的实测发散波谱。
2.如权利要求1所述的利用频率汉克尔变换分离波场的面波勘探方法，其特征在于，所述计算被勘探处在面波震动下的响应结果所对应的实测格林函数，包括：
依据所述面波，得到所述面波中的主动源，所述主动源为人工制造的震源；
依据所述响应结果，得到所述响应结果中的震动数据，所述震动数据为所述被勘探处在所述主动源的震动下产生的震动信息；
将所述震动数据转换到频域，得到频域的所述震动数据；
将所述主动源的震源时间函数转换到频域，得到频域的所述震源时间函数；
依据频域的所述震动数据和频域的所述震源时间函数，计算实测格林函数。
3.如权利要求1所述的利用频率汉克尔变换分离波场的面波勘探方法，其特征在于，所述计算被勘探处在面波震动下的响应结果所对应的实测格林函数，包括：
依据所述面波，得到所述面波中的被动源，所述被动源为非人工制造的震源；
计算所述被勘探处的任意两个监测设备采集到的两个震动数据之间的互相关函数，得到响应结果；
对所述响应结果应用希尔波特变换，得到变换之后的所述响应结果；
将变换之后的所述响应结果作为实部、所述响应结果作为虚部，得到实测格林函数。
4.如权利要求1所述的利用频率汉克尔变换分离波场的面波勘探方法，其特征在于，所述依据所述第一类零阶汉克尔函数、所述第二类零阶汉克尔函数、所述实测格林函数，得到频率矢量波数域格林数谱，包括：
将所述第一类零阶汉克尔函数与所述第二类零阶汉克尔函数之和乘以二分之一，得到贝塞尔函数；
将所述贝塞尔函数、所述实测格林函数与所述距离的乘积关于所述距离做积分，得到频率矢量波数域格林数谱。
5.如权利要求1所述的利用频率汉克尔变换分离波场的面波勘探方法，其特征在于，所述依据所述实测发散波谱，得到针对所述被勘探处所在的实测地层结构，包括：
依据所述实测发散波谱，绘制图像式的实测频散谱；
对所述实测频散谱应用训练后的神经网络模型，得到针对所述被勘探处所在的实测地层结构。
6.如权利要求5所述的利用频率汉克尔变换分离波场的面波勘探方法，其特征在于，所述神经网络模型的训练方式，包括：
计算地层结构模型在面波震动下的响应结果所对应的理论格林函数；
计算所述理论格林函数的核函数；
依据所述核函数的频散特性，绘制理论频散谱；
依据所述地层结构模型和所述理论频散谱，训练所述神经网络模型。
7.一种利用频率汉克尔变换分离波场的面波勘探装置，其特征在于，所述装置包括如下组成部分：
格林函数计算模块，用于计算被勘探处在面波震动下的响应结果所对应的实测格林函数；
波谱分离模块，用于对所述实测格林函数应用频率-汉克尔变换分离出所述实测格林函数所涵盖的实测发散波谱；
地层结构计算模块，用于依据所述实测发散波谱，得到针对所述被勘探处所在的实测地层结构；
所述实测格林函数为所述被勘探处在面波中的主动源或被动源震动下的响应结果所对应的函数，所述对所述实测格林函数应用频率-汉克尔变换分离出所述实测格林函数所涵盖的实测发散波谱，包括：
依据所述实测格林函数，得到与所述实测格林函数所对应的地震波场；
计算所述地震波场的周期与波速之积；
计算所述周期与所述波速之积的倒数与距离的乘积结果，所述距离为所述主动源与所述被勘探处之间的距离，或者，所述距离为相邻监测设备之间的距离，所述监测设备用于监测所述被勘探处的震动数据；
计算所述乘积结果的第一类零阶汉克尔函数和第二类零阶汉克尔函数；
依据所述第一类零阶汉克尔函数、所述第二类零阶汉克尔函数、所述实测格林函数，得到频率矢量波数域格林数谱；
提取所述频率矢量波数域格林数谱所涵盖的实测发散波谱。
8.一种终端设备，其特征在于，所述终端设备包括存储器、处理器及存储在所述存储器中并可在所述处理器上运行的利用频率汉克尔变换分离波场的面波勘探程序，所述处理器执行所述利用频率汉克尔变换分离波场的面波勘探程序时，实现如权利要求1-6任一项所述的利用频率汉克尔变换分离波场的面波勘探方法的步骤。
9.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质上存储有利用频率汉克尔变换分离波场的面波勘探程序，所述利用频率汉克尔变换分离波场的面波勘探程序被处理器执行时，实现如权利要求1-6任一项所述的利用频率汉克尔变换分离波场的面波勘探方法的步骤。
CN202211467007.5A 2022-11-22 2022-11-22 利用频率汉克尔变换分离波场的面波勘探方法及相关装置 Active CN115755176B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211467007.5A CN115755176B (zh) 2022-11-22 2022-11-22 利用频率汉克尔变换分离波场的面波勘探方法及相关装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211467007.5A CN115755176B (zh) 2022-11-22 2022-11-22 利用频率汉克尔变换分离波场的面波勘探方法及相关装置

Publications (2)

Publication Number Publication Date
CN115755176A CN115755176A (zh) 2023-03-07
CN115755176B true CN115755176B (zh) 2023-06-13

Family

ID=85336778

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211467007.5A Active CN115755176B (zh) 2022-11-22 2022-11-22 利用频率汉克尔变换分离波场的面波勘探方法及相关装置

Country Status (1)

Country Link
CN (1) CN115755176B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5588032A (en) * 1992-10-14 1996-12-24 Johnson; Steven A. Apparatus and method for imaging with wavefields using inverse scattering techniques
CN102096100A (zh) * 2010-11-24 2011-06-15 中国石油天然气集团公司 测井曲线与地震记录全波匹配方法及装置
US8448115B1 (en) * 2012-06-18 2013-05-21 Mentor Graphics Corporation Through silicon via impedance extraction

Family Cites Families (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7689949B1 (en) * 2007-03-16 2010-03-30 Cadence Design Systems, Inc. Evaluating Green's functions for modeling multilayer integrated circuits
EP2847574A1 (en) * 2012-05-07 2015-03-18 Helmholtz Zentrum München - Deutsches Forschungszentrum für Gesundheit und Umwelt (GmbH) Apparatus and method for frequency-domain thermo-acoustic tomographic imaging
US9442209B2 (en) * 2012-07-10 2016-09-13 Pgs Geophysical As Methods and systems for reconstruction of low frequency particle velocity wavefields and deghosting of seismic streamer data
CN103593510B (zh) * 2013-10-25 2016-06-29 西安电子科技大学 基于互易性原理的粗糙面与目标复合电磁散射仿真方法
US10074993B2 (en) * 2014-09-11 2018-09-11 Cpg Technologies, Llc Simultaneous transmission and reception of guided surface waves
US10033197B2 (en) * 2015-09-09 2018-07-24 Cpg Technologies, Llc Object identification system and method
US9927477B1 (en) * 2015-09-09 2018-03-27 Cpg Technologies, Llc Object identification system and method
US10498006B2 (en) * 2015-09-10 2019-12-03 Cpg Technologies, Llc Guided surface wave transmissions that illuminate defined regions
CN105184010A (zh) * 2015-09-25 2015-12-23 天津城建大学 基于快速多极间接边界元法的高频地震波散射模拟方法
CA3055749A1 (en) * 2017-03-08 2018-09-13 Saudi Arabian Oil Company Automated system and methods for adaptive robust denoising of large-scale seismic data sets
JP6945895B2 (ja) * 2017-10-12 2021-10-06 サウザン・ユニバーシティ・オブ・サイエンス・アンド・テクノロジー 表面波探査方法および端末デバイス
CN109239712B (zh) * 2018-07-24 2023-03-21 哈尔滨工程大学 基于水下声场和声能流的噪声探测方法
CN111164462B (zh) * 2018-08-06 2022-05-06 南方科技大学 一种人工源面波勘探方法、面波勘探装置及终端设备
CN110208853B (zh) * 2019-05-30 2020-07-14 中国地质大学（北京） 基于自由界面地震波场导数重建的波动方程保幅偏移方法
CN113945968B (zh) * 2021-10-19 2022-04-29 中国矿业大学(北京) 不连续地质体的绕射波成像方法、装置及电子设备
CN114925317B (zh) * 2022-05-17 2023-03-31 北京智芯仿真科技有限公司 一种用于集成电路的快速汉克尔变换方法及装置
CN114722341B (zh) * 2022-05-17 2022-08-12 北京智芯仿真科技有限公司 一种用于集成电路的汉克尔变换滤波器的改进方法及装置
CN115236738A (zh) * 2022-06-27 2022-10-25 中国地质科学院地球物理地球化学勘查研究所 一种基于汉克尔变换的勒夫波频散提取方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5588032A (en) * 1992-10-14 1996-12-24 Johnson; Steven A. Apparatus and method for imaging with wavefields using inverse scattering techniques
CN102096100A (zh) * 2010-11-24 2011-06-15 中国石油天然气集团公司 测井曲线与地震记录全波匹配方法及装置
US8448115B1 (en) * 2012-06-18 2013-05-21 Mentor Graphics Corporation Through silicon via impedance extraction

Also Published As

Publication number Publication date
CN115755176A (zh) 2023-03-07

Similar Documents

Publication Publication Date Title
Foti et al. Guidelines for the good practice of surface wave analysis: a product of the InterPACIFIC project
KR101548976B1 (ko) 지진 표면파들의 파형들을 사용하는 토양 특성들의 추정
EP2030046B1 (en) Vh reservoir mapping
CN104280775B (zh) 一种基于全波形矢量偏移叠加的微地震监测定位方法
CN111538075B (zh) 干热岩勘探方法、装置、电子设备及存储介质
US20120016592A1 (en) Image domain signal to noise estimate with borehole data
WO2012139082A1 (en) Event selection in the image domain
US20120116682A1 (en) Energy density and stress imaging conditions for source localization and characterization
CN112327358B (zh) 一种粘滞性介质中声波地震数据正演模拟方法
Yuan et al. Numerical comparison of time-, frequency-and wavelet-domain methods for coda wave interferometry
CN106094027A (zh) 一种垂直地震剖面vsp钻前压力预测方法和系统
CN111487692A (zh) 一种盐间页岩油韵律层地震响应特征及储层厚度预测方法
CN110850469A (zh) 一种基于克希霍夫积分解的地震槽波深度偏移的成像方法
CN113703044B (zh) 古河道宽度的校正方法、装置、电子设备及存储介质
CN115755176B (zh) 利用频率汉克尔变换分离波场的面波勘探方法及相关装置
Kushnir et al. Determining the microseismic event source parameters from the surface seismic array data with strong correlated noise and complex focal mechanisms of the source
US10317543B2 (en) Estimation of a far field signature in a second direction from a far field signature in a first direction
CN103513279A (zh) 一种基于地震波波动方程的照明分析计算方法及计算装置
Wei et al. A high-resolution microseismic source location method based on contrast source algorithm
US20110276272A1 (en) Time reverse imaging attributes
Hailemikael et al. From ambient vibration data analysis to 1D ground-motion prediction of the Mj 5.9 and the Mj 6.5 Kumamoto earthquakes in the Kumamoto alluvial plain, Japan
Bhaumik et al. Higher-order thin layer method (HTLM) based wavefield modeling approach
CN113552617B (zh) 小尺度缝洞体的量化方法、装置、电子设备及存储介质
Pippin Relative source time functions, spectral ratios, and spall for the source physics experiment phase I explosions
Özkan et al. Coda-derived source properties estimated using local earthquakes in the Sea of Marmara, Türkiye

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