微动中多模式面波频散曲线的映射式提取方法
侯小平
摘 要：折射微动（Refraction Microtremor）面波探勘方法现在已经广泛运用于地震工程勘探中。目前主要通过手动提取来获得该方法的频散曲线，而且手动拾取依赖于处理人员经验，主观性比较强，工作耗时多。本文给出了一种基于最小二乘拟合的频散曲线拾取方法。它主要包含两个步骤：第一步从p-f谱上提取最大梯度点分布；第二步是通过最小二乘拟合得到频散曲线的趋势，完成频散曲线的拾取。
关键词：折射微动；频散曲线；拾取
中图分类号：P631.4 文献标识码：A 文章编号：1671-2064（2018）06-0213-02
1 引言
折射微动（Refraction Microtremor，ReMi）方法是利用背景噪音（比如汽车、机械、人类活动、海浪、风和气候变化等产生的噪音）来探测地下地质结构的一种面波勘探方法[1]。该方法利用多道线性排列进行数据采集，记录时间在30～50s，通常可以探测到的深度为排列长度的1/3～1/2，最大探测深度为排列长度，探测精度与道间距相关。ReMi方法的处理主要分为三个步骤：第一步是将原始数据从空间和时间（x-t）域变换到慢度截距时间（p-τ）域，此时的振幅与慢度和截距时间相关；第二步是在截距时间的方向上做傅里叶变换，将数据从p-τ域变换到p-f域；第三步是在p-f域沿着每一个频率做功率谱比计算。
野外采集的原始地震噪音数据经上述流程处理之后，可以得到类似图1所示的p-f谱。从图中我们可以看到频散曲线的趋势（图1黑色实线）：慢度（相速度）随着频率的增大而增大（减小）；这是因为频散曲线的低频部分对应的是长波长的Rayleigh面波，主要反映了深部地质结构的信息；而高频部分则对应短波长的Raleigh面波，主要反映了浅部的地质结构信息。
根据Louie[1]的研究，微振动数据采集过程中，假设来自各个方向的波的能量都相等，视速度是真速度的100%～125%范围内的能量占40.9%，在视速度大于真速度125%以上的范围内分布着高于59%的能量。如此，在一个特定的频率上，RiMi方法拾取的是有显著能量出现的地方的最小视速度值，整个频散曲线是沿着低速包络进行拾取。根据ReMi方法中的频散曲线拾取步骤，需要在每个频率上需要拾取三个点：首先拾取的是具有低相位速度点，它位于慢度-谱比曲线上谱比值刚好与不相干噪音的低谱比值分离时候的位置；其次拾取最佳猜测速度点，它位于谱比值上升梯度最大的地方；最后拾取的位于谱比值最大时候的位置，以此获得研究场地的真相速度及其变化上下限。Pancha等人[2]通过研究表明，Louie提出的通过拾取沿慢度变化的谱比函数最大梯度点来获得p-f谱中频散曲线的低速包络是合适的。在图1中，频散曲线拾取沿着黑色实线趋势进行拾取。其中，可以看到一个与频散曲线趋势变化相反的能量分布（图1黑色虚线所示），这是由于采样频率低于奈奎斯特（Nyquist）频率而导致的F-K混叠。
2 频散曲线拾取方案
根据折射微动方法频散曲线的拾取原则，我们设计了频散曲线拾取方案，包括如下2项内容：
2.1 提取频散点分布
在p-f谱中，功率谱比值是p和f的函数，我们用F（p，f）来表示功率谱比函数。在频率fi上的功率谱比值可以表示为：
（1）
其中，i表示第i个频率，m是频率离散点数。对于某一频率fi，F（p，fi）只随p的变化而变化，我们可用梯度的变化来描述这种变化的过程，我们在该频率上需拾取的速度点位于F（p，fi）上最大值位置（图2a点）和F（p，fi）的谱比值刚好从不相关噪音谱比值分离时候的位置（图2c点）之间的F（p，fi）上最大梯度点所对应的速度（图2b点），用vst来表示。我们的工作的第一步就是求取每个频率上的vst，然后得到整个p-f谱中的vst分布函数H（f）。理想状态下，H（f）即为我们所需的整条频散曲线，但由实际数据得到的频散曲线H（f）随频率的分布呈锯齿状变化，这种变化在FK混叠区域比较大。在这种情况下，如果直接使用H（f）进行后续的反演工作，将得不到满意的地下横波速度结构，因此需要找到H（f）变化趋势，本文主要通过最小二乘拟合来找到该趋势。
2.2 H（f）最小二乘拟合
在得到H（f）之后，通过（2）式对其进行拟合：
（2）
其中，n为拟合的阶次，n越大，拟合曲线与实际曲线越接近；a为对应的系数。从（2）式可以看到，当我们知道了ai的具体值的时候，我们就知道了h（f）的具体表达式。令最小二乘目标函数I为：
（3）
其中ω为加权系数。（3）式可以转化为求多元函数：
（4）
的极小值点问题，将系数a看作是自变量，通过解方程：
（5）
可获得每一项系数a，然后代入（2）式可以获得h（f）的具体表达式。
3 应用实例
我们将本文所提出的方法运用到某一工区的面波资料。在数据采集过程中，采用24道接收、道间距10m、排列长度为230m、记录长度为50s，时间采样率为2ms，噪音来源主要由车辆、风和人员活动。我们将采集到的数据用ReMi方法处理，得到图3所示的p-f谱，然后再运用本文的方法的第一步，在p-f谱上提取最大梯度点的分布H（f）（图3黑色方框）。从图3中可以看出，H（f）在3-10Hz范围内的分布比较规律；在大于10Hz的区域，受FK混叠的影响H（f）的分布比较乱；根据本文的方法的第二步，我们拾取得到整条频散曲线（图3的灰色曲线），该曲线在大于10Hz的F-K混叠区域能跟随低频段（<10Hz）的趋势。
4 分析讨论
从野外采集的微振動数据经折射微动方法处理之后，频散曲线的趋势会受到FK混叠的干扰，在干扰区域频散趋势不明显；如图3中的p-f图，在2.5-10Hz的区域内，我们能够了看到清晰的频散曲线的趋势；在10-20Hz的区域内，受FK混叠的影响，频散趋势不如2.5-12Hz范围内明显，在这个区域进行人工拾取频散曲线的时候，受到的干扰很大，造成频散曲线拾取的精度的缺失。
在计算最大梯度点的时候，由于低频受FK混叠的干扰小，因此最大梯度点的计算是由低频到高频，根据频散曲线的趋势变化来设计拾取原则。从实际数据中的提取的频散曲线并不是严格的单调变化，通常是呈现锯齿状形态，因此当我们找到某一个频率上的最大梯度点之后，下一个频率上的最大梯度点的范围就可以根据这个点来确定。
在最小二乘拟合的时候，由于低频受混叠影响小，因此在低频的区域建议给较大的权系数，高频区域给较小的权系数，以此来减小FK混叠的影响。关于最小二乘拟合的阶次n，如果n给得过大，拟合曲线就会越接近我们得到的H（f），会受到H（f）局部极值的影响，使得我们得到的频散曲线不够平滑，不能很好的反应出频散曲线的趋势。在兼顾频散曲线的趋势和局部极值的影响之下，我们通过实验发现阶次在7-15比较合适。
5 结语
我们提出了一种针对折射微动方法频散曲线的拾取方法，通过最小二乘法拟合来实现。我们在进行频散曲线拾取的时候，拾取的该频率上的功率谱比-慢度曲线上面的功率谱比值梯度变化最大的时候所对应的速度。
参考文献
[1]Louie J N. Faster， better： shear-wave velocity to 100 meters depth from refraction microtremor arrays[J]. Bulletin of the Seismological Society of America，2001，91（2）：347-364.
[2]Pancha A， Anderson J G， Louie J N， et al. Measurement of shallow shear wave velocities at a rock site using the ReMi technique[J]. Soil Dynamics and Earthquake Engineering， 2008，28（7）：522-535.
