一种有效的适合F-J方法从背景噪声资料中提取多阶频散曲线的台阵划分方法

袁艺,陈晓非

1 北京大学地球与空间科学学院,北京 100871 2 中国地震局第二监测中心,西安 710054 3 广东省地球物理高精度成像技术重点实验室(南方科技大学),深圳 518055 4 南方科技大学地球与空间科学系,深圳 518055

0 引言

可准确、高效提取多阶频散曲线的F-J变换方法(Wang et al.,2019)是基于台阵的方法,研究者们将其应用于世界各个区域来得到高分辨率的地下速度结构(Chen et al.,2022; 李雪燕等,2020; Wu et al.,2020; 吴华礼等,2019; Zhan et al.,2020; Zhang et al.,2022;陈春林等,2023).具体说来,需要先将研究区域的所有台站划分为一个个相互间有重叠的小阵列,每个阵列内所有两两台站的互相关函数的集合(CCFs)经过F-J变换后可得到相应阵列的F-J谱,从中提取出多阶频散曲线进而反演得到对应于该小阵列中心下方的1-D速度模型,最终将所有阵列的1-D速度模型整合得到研究区域的3-D速度模型.使用F-J方法成像的效果很大程度上取决于对台阵的划分方式:一方面,若单个小阵列越小,则整个研究区域的小阵列数目越多,得到的1-D速度模型就越多越密,3-D速度模型的分辨率就越高;另一方面,根据Wang等(2019)的理论分析,对于单个小阵列,其内台间距越小,整个小阵列展布越大,则得到的F-J谱质量越高.若要得到高分辨率成像结果,就需要在这两者间做个权衡.除了Li等(2023)基于Voronoi多边形的分区法,目前还没有其他针对应用F-J方法时如何合理分区的论述,最简单常用的分区方式就是用一个大小形状或台站数目恒定的小阵列沿经纬度依次顺序移动,遍历整个研究区域,这种分区方法适用于台间距比较均匀或台阵下方速度结构变化较平缓的情况;Li等(2023)基于Voronoi多边形的分区法优势在于整个分区过程自动进行,且可定量评估相速度测量的不确定度,但同样适用于地下结构变化比较均匀的区域.若研究区域速度结构复杂多变,台阵间距也疏密不均时,根据具体情况灵活设计适合当地的分区方案更有可能得到高分辨率的结果.本文在实际处理此类地下结构复杂且覆盖台网疏密不均区域数据的过程中摸索出一种使用F-J方法时的独特分区思路,可根据具体情况做调整,具有一定的普适性.

1 方法介绍

对某一研究区域内的台站,经过数据处理的各项常规流程后,我们计算得到了两两台站记录的垂直分量背景噪声的互相关函数集合CCFs.我们首先固定阵列内台站数为30个(可根据研究区域台阵数目及台间距变化),即对区域内每一台站,选取距离其最近的29个台站,检查这30个台站组成的小阵列是否与前面的各小阵列重合,若重合则舍弃,继续对下一个台站重复上述操作;若没有重合,则从互相关集合中提取属于该小阵列的部分,并计算对应的F-J谱I(ω,c).遍历整个研究区域后,会发现不同区域对应的F-J谱差别较大,这是因为从F-J谱中提取的频散曲线反映的就是该处的地下速度结构.我们参考地形,将各小阵列的F-J谱按相似程度分类,从而将整个台阵划分为若干个相互间有重叠的大的子区域,每个子区域内的地下结构变化幅度较小.在每个子区域内,再用更小的固定数目小阵列遍历该子区域.由于各子区域的台站疏密及地下速度结构的差异,需分别确定适合的小阵列内台站数目.我们的原则是,在能得到正常形态F-J谱的前提下尽可能减少小阵列的台站数,以提高最终结果的分辨率.因此对于每一个次级区域,一般需要进行多次尝试才能找到最适合的小阵列数目,在用这样的小阵列遍历子区域进而得到对应的各F-J谱后,提取相应的多阶频散曲线并反演对应于小阵列中心的1-D地下速度模型,最终集合所有的1-D模型得到整个研究区域下方的3-D速度模型(详细分区流程图见图1).

图1 推荐分区方法流程图

2 实例研究

我们下载了IRIS台网从2015年1月1日至2020年12月31日共6年的宽频带垂直分量连续波形数据(BHZ),包括310个台站,跨越整个阿拉斯加及加拿大西部的一小部分,目的是得到整个研究区域下方的三维横波速度结构.阿拉斯加台站分布不均(图2),且地质构造运动强烈.南部位于太平洋板块俯冲区域,且多火山分布,台站布设密集度远大于中部及北部,尤其是北部区域台站分布稀疏.这样的地下速度结构复杂且台阵分布不均的情况非常适合使用我们的分区策略.下面我们分别使用简单的分区方法及推荐的分区方法,对比二者得到的结果.

图2 阿拉斯加及加拿大西北部310个台站分布图

2.1 简单分区方法

首先用一个固定大小的矩形框在研究区域内移动,该矩形框内区域即为小区域.我们根据如下考虑设置矩形框的大小:

(1)阿拉斯加台站分布疏密不均,为了在最稀疏的区域也能得到可靠的F-J谱,须保证框内台站数不能过少,也即矩形框不能太小.

(2)本文研究区域为58°N—72°N,167°W—131°W,纬度跨越14°,经度跨越36°,经度范围远大于纬度.

(3)研究区域纬度高,经度上移动2°基本等同于纬度上移动1°.

综上,本文选择了8°(沿纬度)×10°(沿经度)大小的固定矩形框,从研究区域左上角(纬度最高处)开始,每次向东沿经度移动2°,移到东边界后,向南移动1°并回到西边界,然后继续自西向东移动.重复这个过程,直到矩形框的移动范围已经覆盖整个研究区域.

移动完成后,我们得到了94个小阵列,沿纬度有7行,沿经度最多14列,图3为从其中随机挑选的9个小阵列台站分布图.对每一个小阵列都可通过背景噪声互相关和F-J变换方法得到相应的F-J谱,并可从中提取基阶和高阶频散曲线.图4分别为第3行第5列(0305)及第6行第2列(0602)的小阵列得到的F-J谱(左)及从中提取的多阶频散曲线(右).

图3 8°×10°(即纬度跨度×经度跨度)的矩形框在阿拉斯加顺序移动图示

图4 随机选取的两个小阵列0305和0602对应的F-J谱及从中提取到的多阶频散曲线

2.2 推荐分区方法

按照流程图,我们先用30个台站的小阵列遍历整个研究区域,发现不同区域对应的F-J谱差别很大:东部区域阵列的F-J谱中,基阶频散曲线清晰集中,频带较宽,但多未见高阶频散曲线;中部区域F-J谱中能看到多条高阶频散曲线,但基阶频散曲线较东部区域清晰度下降且频带变窄;在南部区域,F-J谱整体分辨率较低,基阶频散曲线集中度进一步下降,频带更窄,一般不见高阶频散曲线.我们根据这些频散曲线的特征,再结合地形,将整个研究区域的台站分为8个大的次级区域(图5),各个次级区域内的频散曲线有共性.接下来,我们按照同样的方法在各个次级区域内用更小数量的小阵列移动.考虑到各个次级区域台站密度不同,地质特征也有差异,其内更小阵列的台站数最好也分别设置.我们依次在各个次级区域内试验不同的小阵列台站数,数量过少得不到可用的F-J谱,数量过多会导致低分辨率,权衡后确定1区、5区、6区及7区内移动的小阵列固定台站数为15个,其余4个区为18个.图6和图7分别对应1区和3区中随机挑选的小阵列,图(a)为小阵列的台站分布情况,图(b)为小台阵记录数据计算得到的F-J谱,从中可提取相应的频散曲线.一般基阶频散曲线的能量更强,且对反演的影响最大,单独使用基阶反演即可确定下方1-D横波速度的大致形态,而高阶频散曲线则在此基础上对更精细的速度结构作约束.大多数情况下高阶频散曲线能量较基阶弱很多,且有非常多的噪声干扰,为保证提取到的高阶频散曲线的准确性,我们首先从F-J谱中拾取基阶频散点,反演得到一个初步的1-D模型,由该模型计算出的高阶频散曲线来指导真实高阶的提取(见图6及图7右侧F-J谱中的红色点线).所有8个次级区域中,1区提取频散曲线29组,2区提取63组,3区77组,4区54组,5区14组,6区55组,7区56组,8区30组,总共提取378组频散曲线(从同一个F-J谱中提取出的基阶及高阶频散曲线为一组),数量较之前的94组大大增加.

图5 八个次级区域各自的台站分布情况

图6 1区中随机挑选的一个小阵列

图7 3区中随机挑选的一个小阵列

2.3 两种分区方法对比

F-J方法是基于台阵的方法,由一个小阵列得到的F-J谱中提取相应的多阶频散曲线,继而反演出近似对应该小阵列平均中心位置下方的1-D横波速度模型(夏江海,2015;李正波,2020;吴高雄,2020).根据频域上F-J谱I(ω,k)的定义(Wang et al.,2019):

(1)

对实际的任意一个阵列,我们只能通过对有限对台间距r的累加来近似上述积分式(1),所以小阵列展布越大,得到的F-J谱就越准确.另外,对于f-k域上的F-J谱,一方面有(Nyquist,1928; Shannon,1949):

(2)

证明了大展布的阵列能够得到更加准确的F-J谱I(ω,k);另一方面,还可以根据f-k域上的台阵响应函数(ARF)来判断阵列布设的影响(Capon,1969; Rost and Thomas,2002):

(3)

台阵可看作一个空间滤波器,对波数的选择范围只和其几何性质,即阵列展布和台站间距这些因素有关(Nishida et al.,2008; Okada and Suto,2003; Picozzi et al.,2010).阵列展布越大,F-J谱分辨率越高;最小台间距越大,旁瓣越明显.此外,Wang等 (2019)的数值模拟实验发现,对于同样数目的台阵,台站随机排布时得到的F-J谱效果优于线型及三角型.根据以上分析,包括所有台站的阵列展布最大,最小台间距最小,得到的F-J谱最准确,分辨率最高.然而一个阵列只能得到一个F-J谱,对应一条1-D横波速度模型,其可看作是相应大展布阵列下方速度结构的平均,为得到研究区域下方的3-D横波速度模型,我们必须将整个阵列拆分为多个小阵列.小阵列展布越大,单个的F-J谱越准确,但总共得到的小阵列数目越少,且单个1-D速度模型的平均效应强,最终得到的3-D模型分辨率也会较低;而小阵列展布越小,F-J谱准确性下降,但得到的1-D横波速度模型的平均效应也在减弱,模型更接近阵列平均中心点下方的真实结构,因此我们需在单个小阵列F-J谱的准确性和最终3-D横波速度模型的分辨率之间做权衡.简单分区方法单个小阵列展布够大,侧重于各小阵列所得F-J谱的准确性,但同时导致平均效应太强,整体分辨率下降;推荐分区法采用两步分区策略:第一步,将所有台站分为8个大的次级区,各次级区内地下结构差别不大,横向速度变化较小,台站密度也相近;第二步根据每个次级区的具体特征,在每个次级区内再次分区,这不仅考虑到了各次级区内合适的小阵列数目可能会不同,还可以合理地缩小小阵列的展布,得到的F-J谱分辨率虽略有下降,但并不影响提取频散曲线.综上,推荐分区方法的两步分区很适合地下结构复杂且覆盖台阵间距不均的区域,在单个小阵列F-J谱分辨率轻微下降的同时,整个研究区域小阵列的数目大幅上升,最终速度模型的分辨率整体上会有提升.图8是两种分区方法如何影响最终3-D横波速度模型分辨率的简单示意图,图8a中白色三角为推荐分区法中的一个小阵列,黄色原点为该小阵列平均中心位置;图8b画出了所有小阵列的平均中心位置;图8c为简单分区法中的一个小阵列;图8d为简单分区法中所有移动阵列的平均中心位置.每一个黄点都对应一个F-J谱及相应的多阶频散曲线,也即对应一个1-D横波速度模型.图8b相比于图8d,黄点在数量、密度及覆盖范围上都有非常明显的提升,可以预见使用推荐分区方法能提高结果的分辨率.

图8 推荐分区方法(a,b)和简单分区方法(c,d)影响结果分辨率比较图

下面我们对比两种分区方法得到的不同频率(周期)的相速度图.基阶频散曲线对反演的影响最大,可确定所得1-D模型的大致形态,高阶频散曲线是在基阶的基础上对更精细的结构进行约束,所以基阶相速度图的比较已经足够我们确定两种分区方法的优劣.我们首先将研究区域沿经纬度网格化,因为简单分区情况各小阵列平均中心相隔相对较远(图8),网格间距设为0.5°,推荐分区情况网格间距设为0.25°.各网格点某频率上的相速度由该点2°范围内对应频率的所有已知相速度(即图8中各黄色圆点处的相速度值)加权平均得到,权重计算公式如下:

(4)

其中di为待求网格点和第i个已知相速度位置的距离.图9为两种分区方法得到的基阶相速度图及推荐分区方法得到的0.401 Hz处的第一到第三高阶相速度图,将基阶相速度图(图9第一、第二行)与阿拉斯加的地形图(图10)相比较后发现,推荐分区方法中,0.0401 Hz的相速度图上,北部的大范围低速区对应位处阿拉斯加北部的北坡盆地(Colville),盆地南缘对应布鲁克斯岭(Brooks Range),频散图上南部的低速区位于Denali断层以南,处于太平洋板块向北美板块俯冲的俯冲带区域;0.0125 Hz的相速度图上,东北部的高速区对应北美克拉通块体.而这些特征在图9的简单分区基阶相速度图上并没有体现,也很难找到其他任何与阿拉斯加地质地形相对应的点.相速度频散曲线的准确程度直接影响最终横波速度模型的可靠性,显然相较简单分区方法,推荐分区方法得到的频散图与阿拉斯加地质地形情况吻合度更好,分辨率更高.除此以外,图9第三行为推荐分区方法得到对应0.401 Hz的第一到第三高阶相速度图,从其对研究区域的覆盖率可以看出,尽管推荐分区方法减小了单个小阵列的展布,一定程度上降低了其得到的F-J谱的分辨率,但我们依然充分利用了F-J方法的优势,提取到足够的各高阶频散曲线.综上,相较简单分区,推荐分区方法得到的基阶频散曲线更准确,各频率对应的相速度图与阿拉斯加地质特征有更好的对应,可反演得到更准确的地下横波速度模型;推荐分区方法提取的高阶频散曲线至少可达第三高阶,且各高阶相速度图占研究区域的比例不低,能够为已得到的速度模型的精细结构提供约束,提高最终地下横波速度模型的纵向分辨率.

图10 阿拉斯加地形图

3 结论

F-J变换是基于台阵的方法,对一特定的研究区域,不同的台阵划分方案会导致分辨率相差悬殊的结果.本文针对F-J变换提出一种从总台阵到若干个边界相互有重叠的次级区台阵再到各次级区内的更小台阵的二级台阵划分方案,尤其适合于地下结构较复杂且其内台站间距分布不均的区域,充分考虑了研究区域的地域特征,更有可能得到较高分辨率的结果.我们将简单的分区方法和本文推荐的分区方法分别应用于美国阿拉斯加及加拿大西北部区域,在数据完全相同的情况下,推荐分区方法得到的基阶相速度图分辨率远高于简单分区方法的结果,并可与地表明显的地质特征相对应;同时推荐分区充分利用了F-J方法的优势,得到高阶相速度图对研究区域覆盖率较高,可约束地下更为精细的速度结构,提高最终速度模型的纵向分辨率.针对某一具体的研究区域,灵活应用推荐分区方法有望提高最终地下速度模型在纵向及横向的分辨率.