一种基于F-K域的被动源面波成像方法及数据筛选方法

摘要

本发明提供了一种基于F‑K域的被动源面波成像数据筛选方法，包括如下步骤：S1.对原始噪声数据进行预处理，并分割得到不同时间的噪声片段；S2.使用F‑K变换将噪声片段变换到频率波数域，对频率波数域的数据进行归一化，再拟合正负波数域的能量分布；S3.将能量分布与原始能量分布进行相关计算，得到不同时间片段的正负相关系数；S4.对不同时间片段的正负相关系数的大小进行降序排列，得到筛选后的正负因果叠加互相关函数；S5.通过正负因果叠加互相关函数，得到筛选后的数据。本发明可以自动检测具有高信噪比的噪声时间片段，无需人工干预，实现了数据快速自动化地筛选。

权利要求

1.一种基于F‑K域的被动源面波成像数据筛选方法，其特征在于，包括如下步骤：S1.对原始噪声数据进行预处理，并对其进行分割得到不同时间的噪声片段；

S2.使用F‑K变换将所述噪声片段变换到频率波数域，对频率波数域的数据进行归一化，并保留数值大于0.8的能量，再通过最小二乘法拟合正负波数域的能量分布；

S3.将所述能量分布与原始能量分布进行相关计算，得到不同时间片段的正负相关系数；

S4.对所述不同时间片段的正负相关系数的大小进行降序排列，对相关系数数值大于等于0.6的噪声片段进行互相关运算并进行叠加，舍弃相关系数数值小于0.6的片段，得到筛选后的正负因果叠加互相关函数；

S5.通过所述正负因果叠加互相关函数，得到筛选后的数据。

2.根据权利要求1所述的基于F‑K域的被动源面波成像数据筛选方法，其特征在于，在步骤S1中，所述预处理包括去均值、去趋势、重采样、去坏道、带通滤波、时间窗口分段、带通滤波、时域归一化和谱白化处理。

3.根据权利要求2所述的基于F‑K域的被动源面波成像数据筛选方法，其特征在于，所述预处理的过程为：依次进行去均值、去趋势、重采样、去坏道、带通滤波、时间窗口分段、带通滤波、时域归一化和谱白化处理。

4.根据权利要求1所述的基于F‑K域的被动源面波成像数据筛选方法，其特征在于，在步骤S3中，所述正负相关系数的公式为：+ ‑

其中，C和C 分别为正负因果时间段相关系数， 和 分别为归一化F‑K域的负、正波数域数值大于0.8的能量对应的纵坐标， 和 分别为负、正波数域最小二乘拟合直线对应的纵坐标，corr为相关算子。

5.根据权利要求4所述的基于F‑K域的被动源面波成像数据筛选方法，其特征在于，在步骤S4中，所述正负因果叠加互相关函数的公式为：‑ ‑ + +

S (t)为C降序排列前m个时间片段反因果部分的累积叠加互相关函数，S (t)为C降序‑ +排列前n个时间片段因果部分的累积叠加互相关函数，Xi 和Xi 分别为第i个片段互相关函‑ +数的反因果部分和因果部分，取对应的S与S分别作为数据筛选之后叠加互相关函数的反因果部分和因果部分。

6.一种基于F‑K域的被动源面波成像方法，其特征在于，包括：权利要求1～5中任一项所述的基于F‑K域的被动源面波成像数据筛选方法；对S5中筛选后的数据进行频散成像，拾取频散曲线，反演得到近地表横波速度模型。

说明书

一种基于F‑K域的被动源面波成像方法及数据筛选方法

技术领域

[0001] 本发明涉及地质勘测技术领域、被动源面波成像技术领域，尤其是涉及一种基于F‑K域的被动源面波成像方法及数据筛选方法。

背景技术

[0002] 面波分析方法被广泛应用于测量近地表横波波速结构。然而，随着人类对城市地下空间勘探需求的不断增加，且许多城市不宜进行大规模的主动源地震勘探，被动源面波方法得到了迅速发展，特别是在城市近地表调查应用中受到相关从业者的青睐。然而，该方法在实际应用中面临着诸多挑战，比如不均匀的噪声源分布和复杂的场地环境，会严重影响面波的频散成像效果，导致无法拾取准确的频散曲线，导致反演不准或失败。因此，需要进行数据筛选来提高被动源面波频散成像的精度，这一过程通常由人工完成，对于长时间的大型数据集处理来说，既费时又无聊。随着被动源技术在近地表勘探的大规模应用，有学者对近地表被动源面波数据筛选技术进行了研究。目前有基于原始信噪比进行加权叠加、互相关函数不对称性、τ‑p域等被动源数据筛选方法，然而目前的自动筛选技术并不是真正意义上的自动化，它们通常是半自动的，需要在一定的阶段进行人工干预，效率需要进一步提高。

[0003] 由于低成本、无损的特点，被动源面波成像方法在城市近地表勘探中表现出巨大的潜力。在被动源面波成像处理中，分割片段叠加是极为重要的一步。通过将长时间的噪声记录进行分割，然后再对其叠加，可以显著地提高面波能量的成像效果。对于成像质量较好的片段，它在进行叠加时会突显有效信号，而对于成像质量较差的片段，随着它的不断叠加不仅会降低成像效果，而且还会将微弱的有效信号淹没于噪声之中。通过数据筛选对提高被动源面波能量的成像精度具有深远意义。拾取准确的面波频散曲线，并对其进行反演，是近地表横波速度模型建立是否可靠的关键。在近地表被动源面波成像处理中，由于环境噪声的源分布未知和复杂的人为环境，有时很难获得高质量的互相关函数或可靠的面波频散曲线。为了提高面波频散成像的质量，需要进行数据筛选。因此，提供一种基于F‑K域的被动源面波成像方法及数据筛选方法。

发明内容

[0004] 本发明提供一种基于F‑K域的被动源面波成像方法及数据筛选方法，该方法基于噪声数据在F‑K域中面波能量的分布特征，通过最小二乘法拟合能量线性分布，并根据定义的正负因果相关系数对噪声时间片段进行排序，自动检测具有高信噪比的噪声时间片段，无需人工干预。

[0005] 本发明实施例的一方面公开了一种基于F‑K域的被动源面波成像数据筛选方法，包括如下步骤：

[0006] S1.对原始噪声数据进行预处理，并对其进行分割得到不同时间的噪声片段；

[0007] S2.使用F‑K变换将所述噪声片段变换到频率波数域，对频率波数域的数据进行归一化，并保留数值大于0.8的能量，再通过最小二乘法拟合正负波数域的能量分布；

[0008] S3.将所述能量分布与原始能量分布进行相关计算，得到不同时间片段的正负相关系数；

[0009] S4.对所述不同时间片段的正负相关系数的大小进行降序排列，对相关系数数值大于等于0.6的噪声片段进行互相关运算并进行叠加，舍弃相关系数数值小于0.6的片段，得到筛选后的正负因果叠加互相关函数；

[0010] S5.通过所述正负因果叠加互相关函数，得到筛选后的数据。

[0011] 在一些实施例中，在步骤S1中，所述预处理包括去均值、去趋势、重采样、去坏道、带通滤波、时间窗口分段、带通滤波、时域归一化和谱白化处理。

[0012] 在一些实施例中，所述预处理的过程为：依次进行去均值、去趋势、重采样、去坏道、带通滤波、时间窗口分段、带通滤波、时域归一化和谱白化处理。

[0013] 在一些实施例中，在步骤S3中，所述正负相关系数的公式为：

[0014]+ ‑

[0015] 其中，C和C分别为正负因果时间段相关系数， 和 分别为归一化F‑K域的负、正波数域数值大于0.8的能量对应的纵坐标， 和 分别为负、正波数域最小二乘拟合直线对应的纵坐标，corr为相关算子。

[0016] 在一些实施例中，在步骤S4中，所述正负因果叠加互相关函数的公式为：

[0017]

[0018] S‑(t)为C‑降序排列前m个时间片段反因果部分的累积叠加互相关函数，S+(t)为C+‑ +降序排列前n个时间片段因果部分的累积叠加互相关函数，Xi和Xi分别为第i个片段互相‑ +
关函数的反因果部分和因果部分，取对应的S与S分别作为数据筛选之后叠加互相关函数的反因果部分和因果部分。

[0019] 本发明实施例的另一方面公开了一种基于F‑K域的被动源面波成像方法，包括：

[0020] 上述中任一项所述的基于F‑K域的被动源面波成像数据筛选方法；

[0021] 对S5中筛选后的数据进行频散成像，拾取频散曲线，反演得到近地表横波速度模型。

[0022] 综上所述，本发明至少具有以下有益效果：

[0023] 本发明基于噪声数据在F‑K域中面波能量的分布特征，通过最小二乘法拟合能量线性分布，并根据定义的正负因果相关系数对噪声时间片段进行排序，自动检测具有高信噪比的噪声时间片段，无需人工干预，实现了数据快速自动化地筛选，并显著地提高了被动源面波频散成像的精度，有效拓宽了面波能量频带。

附图说明

[0024] 为了更清楚地说明本发明实施例技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

[0025] 图1为本发明中所涉及的基于F‑K域的被动源面波成像数据筛选方法的步骤示意图。

[0026] 图2为本发明中所涉及的基于F‑K域的被动源面波成像方法的流程示意图。

[0027] 图3a为本发明中所涉及的A地区的数据采集示意图。

[0028] 图3b为本发明中所涉及的A地区的采集时间记录的示意图。

[0029] 图4a～图4c为本发明中所涉及的A地区的正因果质量好片段F‑K谱、F‑K域能量拟合效果和成像结果的示意图。

[0030] 图4d～图4f为本发明中所涉及的A地区的正因果质量差片段F‑K谱、F‑K域能量拟合效果和成像结果的示意图。

[0031] 图5a～图5c为本发明中所涉及的A地区的负因果质量好片段F‑K谱、F‑K域能量拟合效果和成像结果的示意图。

[0032] 图5d～图5f为本发明中所涉及的A地区的负因果质量差片段F‑K谱、F‑K域能量拟合效果和成像结果的示意图。

[0033] 图6a～图6c为本发明中所涉及的A地区的正因果、负因果和正负因果筛选前成像结果的示意图。

[0034] 图6d～图6f为本发明中所涉及的A地区的正因果、负因果和正负因果筛选后成像结果的示意图。

[0035] 图7a为本发明中所涉及的B地区的数据采集示意图。

[0036] 图7b为本发明中所涉及的B地区的采集时间记录的示意图。

[0037] 图8a～图8c为本发明中所涉及的B地区的正因果、负因果、正负因果筛选前成像结果的示意图。

[0038] 图8d～图8f为本发明中所涉及的B地区的正因果、负因果、正负因果筛选后成像结果的示意图。

[0039] 图9a为本发明中所涉及的B地区的4‑20Hz主动源面波频散成像结果的示意图。

[0040] 图9b为本发明中所涉及的B地区的4‑60Hz主动源面波频散成像结果的示意图。

[0041] 图9c为本发明中所涉及的B地区的主被动联合对比的示意图。

具体实施方式

[0042] 在下文中，仅简单地描述了某些示例性实施例。正如本领域技术人员可认识到的那样，在不脱离本发明实施例的精神或范围的情况下，可通过各种不同方式修改所描述的实施例。因此，附图和描述被认为本质上是示例性的而非限制性的。

[0043] 下文的公开提供了许多不同的实施方式或例子用来实现本发明实施例的不同结构。为了简化本发明实施例的公开，下文中对特定例子的部件和设置进行描述。当然，它们仅仅为示例，并且目的不在于限制本发明实施例。此外，本发明实施例可以在不同例子中重复参考数字和/或参考字母，这种重复是为了简化和清楚的目的，其本身不指示所讨论各种实施方式和/或设置之间的关系。

[0044] 下面结合附图对本发明的实施例进行详细说明。

[0045] 如图1所示，本发明实施例的一方面公开了一种基于F‑K域的被动源面波成像数据筛选方法，包括如下步骤：

[0046] S1.对原始噪声数据进行预处理，并对其进行分割得到不同时间的噪声片段。

[0047] 在一些实施例中，在步骤S1中，所述预处理包括去均值、去趋势、重采样、去坏道、带通滤波、时间窗口分段、带通滤波、时域归一化和谱白化处理。

[0048] 在一些实施例中，所述预处理的过程为：依次进行去均值、去趋势、重采样、去坏道、带通滤波、时间窗口分段、带通滤波、时域归一化和谱白化处理。

[0049] S2.使用F‑K变换将所述噪声片段变换到频率波数域，对频率波数域的数据进行归一化，并保留数值大于0.8的能量，再通过最小二乘法拟合正负波数域的能量分布。

[0050] S3.将所述能量分布与原始能量分布进行相关计算，得到不同时间片段的正负相关系数。

[0051] 在一些实施例中，在步骤S3中，所述正负相关系数的公式为：

[0052]

[0053] 其中，C+和C‑分别为正负因果时间段相关系数， 和 分别为归一化F‑K域的负、正波数域数值大于0.8的能量对应的纵坐标， 和 分别为负、正波数域最小二乘拟合直线对应的纵坐标，corr为相关算子。

[0054] S4.对所述不同时间片段的正负相关系数的大小进行降序排列，对相关系数数值大于等于0.6的噪声片段进行互相关运算并进行叠加，舍弃相关系数数值小于0.6的片段，得到筛选后的正负因果叠加互相关函数。

[0055] 在一些实施例中，在步骤S4中，所述正负因果叠加互相关函数的公式为：

[0056]

[0057] S‑(t)为C‑降序排列前m个时间片段反因果部分的累积叠加互相关函数，S+(t)为C+‑ +降序排列前n个时间片段因果部分的累积叠加互相关函数，Xi和Xi分别为第i个片段互相‑ +
关函数的反因果部分和因果部分，取对应的S与S分别作为数据筛选之后叠加互相关函数的反因果部分和因果部分。

[0058] S5.通过所述正负因果叠加互相关函数，得到筛选后的数据。

[0059] 如图1和图2所示，本发明实施例的另一方面公开了一种基于F‑K域的被动源面波成像方法，包括：

[0060] 上述中任一项所述的基于F‑K域的被动源面波成像数据筛选方法；

[0061] 对S5中筛选后的数据进行频散成像，拾取频散曲线，反演得到近地表横波速度模型。

[0062] 综上所述，本发明的发明构思如下：

[0063] 对野外噪声数据进行预处理之后进行分段，将其变换到F‑K域。由于面波频散能量对应在F‑K域表现为一条倾斜的直线，通过对正负波数域的能量进行最小二乘线性拟合得到能量的主要分布趋势，再计算原始能量与拟合直线的正负相关系数。当F‑K域面波能量越连续，占比能量越多时，拟合的直线与F‑K域能量的分布趋势越接近，其相关系数就越高。当F‑K域能量杂乱无章，面波能量断断续续，拟合的直线对F‑K域能量的分布趋势表达就越低，其相关系数就越低，这样通过相关系数的高低就可以完成高低信噪比片段的分类。为提高相关系数计算效率，对F‑K域数据进行归一化，将数值大于0.8的能量当作有效能量参与计算。然后将时间片段分别按照正负相关系数进行排序，得到筛选后叠加互相关函数，最后进行频散成像。至此就完成了整套数据自动筛选和频散成像过程，该方法操作简单，无需人为经验给出阈值和过多参数。

[0064] 上述过程简述如下：

[0065] 1)对野外噪声数据进行一系列预处理，对其进行分割得到多个时间片段的噪声数据；

[0066] 2)将时间空间域噪声数据片段进行F‑K变换，变换到频率波数域，对频率波数域的数据进行归一化，并保留数值大于0.8的能量，最后通过最小二乘法拟合正负波数域的能量分布；

[0067] 3)将拟合的能量分布与原始能量分布进行相关运算，得到所有片段的正相关系数C+和负相关系数C‑；

[0068] 4)将所有噪声片段按照C+和C‑的大小进行降序排列，对C+和C‑较大的噪声片段进行互相关运算并进行叠加，舍弃C+和C‑较小的片段，最后得到筛选后的正负因果叠加互相关函数。

[0069] 通过以上，实现了利用噪声数据在F‑K域的不同能量特征进行被动源面波成像数据自动筛选。

[0070] 为了进一步说明本发明的技术方案，举例如下：

[0071] 以某西部A地区城市和B地区近地表被动源资料的筛选过程和反演结果为例，进行分析。

[0072] 选择了毗邻城市公路和建筑施工的2H时长背景噪声数据作为典型实例进行分析。A地区采集示意图如图3(a)所示，其原始记录波形如图3(b)所示,在本例中，使用了31个垂直分量的检波器数据，选取了10S作为合适的时间分割长度，将第一道检波器作为虚震源，并采用了干涉法进行虚炮集的恢复。然后利用图2的流程对不同实际噪声片段进行分类得到如图4a～4f和图5a～5f不同质量的噪声片段，最后将相关系数高的片段进行叠加得到筛选后的频散成像结果。

[0073] 通过图6a～6f的对比不难发现，在完成数据筛选之后，无论是反因果部分还是因果部分的能量成像，都比数据筛选前的成像效果要好。在筛选前可接受的面波有效频带为5‑14Hz,而在筛选后可接受的面波有效频带为5‑20Hz。不仅如此，而且对于筛选之后的面波能量具有更好的分辨能力和连续性，高阶能量得以凸显，有利于后续反演多解性的降低和反演精度的提升，说明了方法的有效性。

[0074] 为了进一步说明方法效果，对B地区数据进行筛选。B地区采集示意图如图7(a)所示。使用了31个垂直分量的检波器的采集数据，道间距为5m，其波形如图7(b)所示。从中可以看出多次出现强振幅噪声，主要噪声源为路过的汽车，因此该数据是典型的交通高频背景噪声。对图7(b)的数据进行一系列预处理，以第一道数据为虚拟源，在完成预处理之后，便开始按照图2的自动筛选流程完成数据筛选，得到如图8a～8f的筛选结果。可以看出经过数据筛选后，扩宽了面波的有效频带，改善成像结果，有利于后续频散曲线提取和反演。为了说明筛选数据的可靠性，在相同区域进行主动源采集,通过图9a～9c的对比，主被动源的频散能量吻合较好，同时筛选后的成像结果补充了主动源的低频带能量，具有较好的保真性。

[0075] 通过两个不同地区的数据实验发现，使用基于F‑K域能量特征的被动源面波成像数据自动筛选方法，拓宽了面波能量频带，可以改善面波频散成像的质量获得较为准确的面波频散曲线。经过与主动源资料对比发现，该方法具有较好的准确性和实用性，为后续的反演提供了可靠的数据。

[0076] 综上，本发明的主要在于：

[0077] 方法原理的可靠性。主要使用理论基础较为成熟、较为简单的F‑K变换、最小二乘法、相关等方法。充分利用面波频散能量对应在F‑K域表现为能量聚焦的倾斜直线的特征，通过最小二乘法拟合F‑K域正负波数域的能量线性分布，计算与原始能量的相关程度来表征F‑K域能量的离散程度和数据质量的好坏，在原理上具有较高的可靠性。同时，相比其他数据变换法，F‑K变换在任何情况下都能较好地展现面波的能量特征。

[0078] 操作简单易实现。F‑K变换、最小二乘法、相关等算法代码容易获取，通过最小二乘法的拟合结果与原始能量相关得到正负相关系数，很容易就能将不同质量的噪声片段进行分类。直接保留相关系数较大的片段并对其进行叠加成像，即可实现被动源面波成像数据自动筛选。

[0079] 筛选结果的可靠性。在筛选过程中，首先对原始噪声数据进行了一系列预处理尽可能保留了所关心的面波数据，充分考虑面波能量在F‑K域的聚焦特征，保证了数据分类与筛选的可靠性，不但适用于检波器采集的数据，也适用于DAS采集的数据。相比较其他面波反演方法，该方法计算参数少、无需太多人工干预，是一种简单、高效的数据筛选方法。

[0080] 以上所述实施例是用以说明本发明，并非用以限制本发明，所以举例数值的变更或等效元件的置换仍应隶属本发明的范畴。

[0081] 由以上详细说明，可使本领域普通技术人员明了本发明的确可达成前述目的，实已符合专利法的规定。

[0082] 尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，应当指出的是，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

[0083] 应当注意的是，上述有关流程的描述仅仅是为了示例和说明，而不限定本说明书的适用范围。对于本领域技术人员来说，在本说明书的指导下可以对流程进行各种修正和改变。然而，这些修正和改变仍在本说明书的范围之内。

[0084] 上文已对基本概念做了描述，显然，对于阅读此申请后的本领域的普通技术人员来说，上述发明披露仅作为示例，并不构成对本申请的限制。虽然此处并未明确说明，但本领域的普通技术人员可能会对本申请进行各种修改、改进和修正。该类修改、改进和修正在本申请中被建议，所以该类修改、改进、修正仍属于本申请示范实施例的精神和范围。

[0085] 同时，本申请使用了特定词语来描述本申请的实施例。例如“一个实施例”、“一实施例”、和/或“一些实施例”意指与本申请至少一个实施例有关的某一特征、结构或特性。因此，应当强调并注意的是，本说明书中在不同位置两次或以上提及的“一实施例”或“一个实施例”或“一替代性实施例”并不一定是指同一实施例。此外，本申请的一个或多个实施例中的某些特征、结构或特点可以进行适当的组合。

[0086] 此外，本领域的普通技术人员可以理解，本申请的各方面可以通过若干具有可专利性的种类或情况进行说明和描述，包括任何新的和有用的过程、机器、产品或物质的组合，或对其任何新的和有用的改进。因此，本申请的各个方面可以完全由硬件实施、可以完全由软件(包括固件、常驻软件、微代码等)实施、也可以由硬件和软件组合实施。以上硬件或软件均可被称为“单元”、“模块”或“系统”。此外，本申请的各方面可以采取体现在一个或多个计算机可读介质中的计算机程序产品的形式，其中计算机可读程序代码包含在其中。

[0087] 本申请各部分操作所需的计算机程序代码可以用任意一种或以上程序设计语言编写，包括如Java、Sca l a、Sma l l ta l k、Eiffe l、JADE、Emera l d、C++、C#、VB.NET、Python等的面向对象程序设计语言、如C程序设计语言、Vi sua l Bas ic、Fortran2103、Per l、COBOL2102、PHP、ABAP的常规程序化程序设计语言、如Python、Ruby和Groovy的动态程序设计语言或其它程序设计语言等。该程序代码可以完全在用户计算机上运行、或作为独立的软件包在用户计算机上运行、或部分在用户计算机上运行部分在远程计算机运行、或完全在远程计算机或服务器上运行。在后种情况下，远程计算机可以通过任何网络形式与用户计算机连接，比如局域网(LAN)或广域网(WAN)，或连接至外部计算机(例如通过因特网)，或在云计算环境中，或作为服务使用如软件即服务(SaaS)。

[0088] 此外，除非权利要求中明确说明，本申请所述处理元素和序列的顺序、数字字母的使用、或其他名称的使用，并非用于限定本申请流程和方法的顺序。尽管上述披露中通过各种示例讨论了一些目前认为有用的发明实施例，但应当理解的是，该类细节仅起到说明的目的，附加的权利要求并不仅限于披露的实施例，相反，权利要求旨在覆盖所有符合本申请实施例实质和范围的修正和等价组合。例如，尽管上述各种组件的实现可以体现在硬件设备中，但是它也可以实现为纯软件解决方案，例如，在现有服务器或移动设备上的安装。

[0089] 同理，应当注意的是，为了简化本申请披露的表述，从而帮助对一个或多个发明实施例的理解，前文对本申请的实施例的描述中，有时会将多种特征归并至一个实施例、附图或对其的描述中。然而，本申请的该方法不应被解释为反映所申明的客体需要比每个权利要求中明确记载的更多特征的意图。相反，发明的主体应具备比上述单一实施例更少的特征。