鄂尔多斯及邻区基于程函方程的面波层析成像

蔡光耀， 王未来*， 吴建平,2， 房立华,2

1 中国地震局地球物理研究所， 北京 100081 2 中国地震局地震观测与地球物理成像重点实验室， 北京 100081

0 引言

鄂尔多斯及其邻区位于华北克拉通的西部，其中鄂尔多斯盆地具有稳定的克拉通基底的特性.鄂尔多斯地块受到来自青藏高原北东向的挤压作用，块体的四周被一系列的断陷盆地带所围，其地块内部比较稳定，而周围盆地带的运动比较复杂(张培震等，2002).鄂尔多斯块体的周围分布着银川—吉兰泰断陷带、河套断陷带、山西断陷带以及渭河断陷带，断陷带一般由若干条断裂带组成，这些断陷带的形成时期和运动方式不完全相同(国家地震局，1988；邓起东等，1999).鄂尔多斯盆地的周围分布着强震活动带，其盆地周缘发生过多次大地震，例如1303年的洪洞大地震、1556年的华县大地震、1695年的临汾大地震、1739年的平罗大地震、1920年的海原大地震等等(顾功叙，1983).

鄂尔多斯地块因其独特的地质条件，引起较多的地质学家和地球物理学家对该区域的研究.大尺度体波走时及远震体波走时层析成像揭示，鄂尔多斯下方的高速异常可以延续至200～300 km深度，大同和山西地堑之下的低速异常可一直延伸至400～500 km(Huang and Zhao，2006；Sun et al.，2008；Zhao et al.，2009；Huang et al.，2009；Lei，2012；Tian and Zhao，2013；Wei et al.，2016；郭慧丽等，2017；陈兆辉等，2018).区域地震波走时层析成像揭示，青藏高原东北缘和鄂尔多斯块体西缘地区从地壳到上地幔顶部范围内的速度结构存在显著的横向不均匀性(杨婷等，2012；王椿镛等，2015；Cheng et al.，2016；Sun and Kennett，2016；汪晟等，2017；郭慧丽和丁志峰，2018).地震面波和噪声层析成像(An and Shi，2006；Zheng et al.，2008,2011；何正勤等，2009；房立华等，2009；黄忠贤，2011；Wang et al.，2012；郑现等，2012；Tang et al.，2013；Legendre et al.，2014；Bao et al.，2015；范文渊等，2015；Shen et al.，2016)揭示了鄂尔多斯中下地壳及上地幔岩石圈具有较快的 S 波速度.由于受地震台站分布不均匀和密集程度的影响，有关鄂尔多斯块体内部及其周边拉张盆地的速度横向变化仍然只能获得粗略的图像，很难得到更加精细的深部结构的变化特征.

中国地震局利用“地震科学探测台阵系统”的观测设备在鄂尔多斯及邻区进行了地震台阵观测.目前已有较多科研人员利用这些地震台阵观测的数据，运用不同的方法对鄂尔多斯及邻区的深部结构进行研究(郭慧丽等，2017；潘佳铁等，2017；Wang et al.，2017；王兴臣等，2017；钟世军等，2017；高翔等，2018).近年来，与鄂尔多斯相邻的华北地区增设了大量的流动台站，台站的覆盖性更好，分布更广泛，数据更加丰富.本文充分利用固定台站和流动台站的观测资料，使用基于程函方程面波成像的方法(Lin et al.，2009；Jin and Gaherty，2015；钟世军等，2017)来研究鄂尔多斯及邻区的深部结构特征，进一步探讨鄂尔多斯盆地南北部岩石圈的差异以及大同火山活动相关的动力学问题.

1 数据

图1 研究区台站分布及区域构造特征红色实线代表构造边界.三角形是宽频带台站，正方形是甚宽频带台站.红色符号代表中国地震台阵探测二期246个流动台站；蓝色符号代表中国地震台阵探测三期309个流动台站；黑色符号代表中国地震台网131个固定台站； 绿色符号表其它来源162个流动台站.Fig.1 Distribution of stations and regional tectonic features in the study areaRed full lines represent tectonic boundaries. Triangles are broadband stations and squares are very broadband stations.Red symbols represent the 246 temporary stations from small arrays of the ChinArray II， blue symbols represent the 309 temporary stations from small arrays of the ChinArray III， black symbols represent the 131 permanent stations of the China Seismic Network， green symbols represent the 162 temporary stations from other sources.

图2 地震事件分布图红色三角形为研究区中心；不同颜色的实心圆圈代表挑选出的不同震级的地震事件.Fig.2 The distribution of seismic eventsThe red pentacle denotes center of research region， solid circles of different colors represent selected seismic events of different magnitudes.

本文使用的数据是“中国地震台阵探测”二期和三期以及其他布设的流动台站再加上固定台站记录的地震数据.本文总共使用了848个地震台站，其中固定台站有131个，流动台站则包括“中国地震台阵探测”二期246个和“中国地震台阵探测”三期309个及其他流动观测台站162个.研究区内的流动地震仪器主要包括两类：宽频带地震仪(地震计CMG-3ESPC和数采Reftek-130，频带范围60 s～50 Hz)和甚宽频带地震仪(地震计CMG-3T和数采Reftek-130 S，频带范围120 s～50 Hz).流动台站和固定台站的分布位置较好的覆盖整个鄂尔多斯地块及其邻边地区(图1).迄今为止研究区的台站分布最完整、最密集，台站间距约35 km.首先得到2006年10月至2018年5月之间的地震事件的波形数据，然后将满足震中距在10°～160°的范围内、震源深度小于50 km，震级大于MS5.5条件下的事件波形数据筛选出来.图2给出了经过筛选得到的657个地震事件相对于台阵中心的位置分布.将筛选得到的地震波形数据进行重采样(1 Hz)、去仪器响应、去均值、去倾斜分量、带通滤波等处理.在处理中，采用了3个频段频率范围分别是：12～20 s、25～55 s、60～150 s，带通滤波参数分别为8～25 s、20～60 s、60～200 s，分别参考相速度为3.4 km·s-1、4 km·s-1、4.3 km·s-1，最大台站对间距分别为75 km、150 km和150 km，光滑长度分别为0.25°.数据的初步处理可以减小一系列的误差，有利于得到分辨率较高的成像结果.

2 方法

Lin等(2009)提出了一种高频近似条件下程函方程的面波层析成像方法，并且使用美国西部地区的噪声数据得到了分辨率较高的相速度分布图.基于程函方程面波层析成像的新方法不经过反演，直接测量不同频率下面波的走时场，通过计算走时场的梯度，然后确定不同周期空间慢度场的分布，最后得到某一频率下面波相速度的二维分布(冯吉坤等，2014).Jin和Gaherty(2015)提出了利用密集台阵阵列反演面波的相速度，在程函方程面波层析成像的基础上进行改进并且开发了基于多通道互相关的面波相速度自动测量系统(ASWMS)，通过该系统使用2006年到2014年美国台阵记录到的800多个事件的宽频带地震数据，得到美国大陆20～100 s之间的瑞雷波的相速度分布图.其基本流程可分为两大部分：台站间的相延迟测量；通过程函方程进行相速度的计算.

2.1 台站间相延迟

Gee和Jordan(1992)提出了广义地震数据泛函数，随后被广泛应用于地壳和上地幔模型的构建.广义地震数据泛函数被应用到多通道的分析上，使用实际观测到的波形记录数据进行互相关计算来提取参数.在时间域内选择需要的波形记录的范围，S1与WSS2的互相关定义为

C(t)=S1*WSS2,

(1)

C(t)表示两相邻台站对的互相关函数；两相邻台站对的波形记录分别为：S1和S2；WS表示从垂向地震记录中隔离出基阶瑞利面波窗函数；WSS2表示应用窗函数WS分离出波形记录S2中面波的主要能量和大部分的尾波，文中图3(带通滤波参数是20～60 s；地震信息：2018-03-26 09∶51∶00.7, -5.462,-151.396,40,MW6.6)给出了一个典型的隔离窗函数WS的位置.采用频时分析的方法估算每个台的所有频段的群延迟时间，以群延迟时间前的2个周期为窗口的起始时间，向后5个周期为窗口的终止时间.

图3 隔离窗口函数WS的位置两条红线表示窗口的起止位置.带通滤波参数是20～60 s，地震信息：2018-03-26 09∶51∶00.7, -5.462, -151.396, 40, MW6.6.Fig.3 Isolate the location of the window function WSThe two red lines represent the starting and ending positions of the window functions.The band pass filter parameter is 20～60 s, seismic information :2018-03-26 09∶51∶00.7, -5.462, -151.396, 40, MW6.6.

五参数高斯小波函数(Gee and Jordan，1992)可以很好地将经过高斯窄带滤波后的互相关函数拟合出来，其中高斯小波函数等于高斯包络函数和余弦函数的乘积，通常表示为

Fi*WcC(t)≈AGa[σ(t-tg)]cos[ω(t-tp)],

(2)

Fi表示第i个中心频率对应的高斯窄带滤波器；Wc表示作用于互相关函数C(t)的窗函数；A为幅度因子；Ga为高斯函数；σ为半个频宽；ω为窄带波形的中心频率；tg和tp分别表示群延迟和相延迟时间.

船舶尺度对船舶领域有很大的影响，尤其是船舶长度。为研究船舶尺度对船舶领域的影响，选取研究水域内除渡船以外的70～90 m和90～110 m两类船舶的船舶领域进行对比分析。同时，船舶在航道中航行时分为靠近沙市的下行方向和靠近埠河的上行方向，通常上行船舶航速小于下行船舶航速，为研究航速与船舶领域之间的关系，将每一种尺寸的船舶分为上行船舶、下行船舶。

2.2 相速度计算

通过慢度向量反演计算相速度的方法，可以压制高次波、低幅度的相速度变化，有利于提高最终相速度图像的分辨率(Jin and Gaherty，2015).相邻两台站间的相延迟时间δτp和慢度矢量之间的关系可以用矢量场积分表示为

(3)

S(r)是慢度矢量；ri表示第i个相邻台站之间的球面路径.将(3)式写成离散的形式：

(4)

SRi和STi分别表示慢度矢量在i处的径向分量和切向分量；drRi和drTi分别表示相邻台站之间的球面路径离散化第i段处的径向分量和切向分量.将反演慢度向量的目标函数表示为

(5)

式中等号右边第一项是观测相延迟与预测相延迟之间的失配值；第二项是慢度光滑的约束，λ表示光滑因子，其中λ和频率有关，本文各周期的光滑因子分别为：3*[0.5，0.4，0.3，0.2，0.2，0.2，0.2，0.2，0.2].反演中各个地震的波前面矢量是分别计算的，以此获得相应格点的相速度分布，然后对各地震的格点相速度结果进行叠加，本文只平均叠加分析平均各方位的平均值，没有考虑或者分析各向异性的分布情况.在叠加过程中，三个周期范围的可容忍相速度范围分别为2.5～4 km·s-1、3～4.5 km·s-1、3.5～5 km·s-1.在图4中给出了4个典型周期的程函反演中用到的台站对，显示了各台站的使用情况和台站对射线方位覆盖情况.

3 结果

本文利用基于程函方程的面波层析成像的新方法，得到鄂尔多斯及邻区12～150 s周期的相速度，图5中展示了12 s、16 s、20 s、30 s、45 s、60 s、90 s、120 s、150 s周期的相速度图像，在一定程度上很好地揭示了深部结构的横向变化特征.基阶瑞利波的相速度大约对1/3波长深度的S波比较敏感，我们计算了不同周期、不同深度的频散灵敏度系数曲线(图6)，其中一维初始速度模型(图6a)参考了Crust1.0.根据频散灵敏度系数曲线易知，不同周期的相速度图像很好的与地表沉积层、地壳、上地幔等不同深部结构的特征对应起来.

图4 台间相延迟测量密度分布图红色三角形代表台站位置，黑色实线代表程函反演中用到的台站对.Fig.4 Density distribution diagram of interstage phase delay measurementThe red triangles represent the station locations and the black solid lines represent the station-pair links used in the equation inversion.

图5 不同周期的相速度成像结果Fig.5 Results of phase velocity imaging in different periods

图6 (a) 计算频散灵敏度系数所用的参考模型； (b)、(c)、(d) 分别表示不同周期范围的基阶Rayleigh面波相速度对横波速度结构的频散灵敏度系数Fig.6 (a) Reference model used in dispersion sensitivity coefficient calculation; (b)、(c)、(d) Respectively represent the dispersion sensitivity coefficients of the phase velocities of fundamental-order Rayleigh surface waves in different period ranges to the shear wave velocity structures

图7 不同周期的相速度成像结果误差分布图Fig.7 Error distribution of phase velocity imaging results in different periods

为了更直观的显示相速度结果的可靠性，我们画了相速度的误差统计图(图7).在12 s的图像中，河套盆地附近存在较大的误差.短周期面波资料不够发育信噪比相对较差；沉积层区域的台站噪声相对较高；河套盆地位于研究区边缘地区，台站对覆盖不好，因此，该地区更容易受到影响.整体来讲，12～60 s相速度的误差比较小能够很好的反映地壳和上地幔深部精细结构的特征，随着周期的增加，由于受到波长的影响导致相速度的误差会有所增加，但是在同一周期研究区域内的相速度误差呈均匀分布，这说明在一定程度上仍能很好的反映地下深部结构的特征.

根据频散灵敏度系数曲线可知，周期为12～20 s的图像主要对深度在12～34 km范围的S波速度敏感，能够反映地表沉积层和地壳的深部结构特征.河套盆地和银川盆地呈现出明显的低速异常，和钟世军等(2017)得到的结果比较吻合.鄂尔多斯盆地东南部的相速度高于西北部，这和郑现等(2012)利用背景噪声得到的结果是一致的.随着周期的增加，由于各盆地沉积层厚度的不同导致低速异常的变化幅度不同.山西断陷盆地由于具有大量的沉积层而显示低速；两侧的太行山和吕梁山具有较浅或者出露的基岩而显示较高的速度分布(周铭等，2016；王霞等，2019).短周期相速度的高低速异常的分布特征很好的将地形地貌的特征刻画出来.

周期为30～60 s的图像主要对深度在28～136 km范围内的S波速度敏感，能够反映下地壳和上地幔的深部结构特征.大同区域的中下地壳存在大面积低速异常的现象，这和郭震等(2015)利用噪声面波和布格重力异常联合反演得到的结果是一致的.大同火山区的低速异常从中下地壳到上地幔的顶部进行不断的延伸和扩展，低速异常变得越来越明显且范围不断扩大，与已有的研究结果一致(Wang et al.，2012；杨婷等，2012；宋美琴等，2013；王霞等，2019).利用接收函数得到鄂尔多斯块体的地壳厚度变化比较平缓约为44 km，华北克拉通中部地壳厚度相对鄂尔多斯地块较薄，在鄂尔多斯块体及其周围整体上表现出西部厚，东部薄，北部厚，南部薄的分布特点(Wang et al.，2017；娄辛辉等，2017).根据成像结果可知，相速度的图像在一定程度上大致反映出研究区地壳厚度差异的分布特征.山西断陷带以北纬38°为分界线在下地壳和上地幔的南北部存在明显的速度差异的分布特征，与现有结果具有较好的一致性(宋美琴等，2013；周铭等，2016；王霞等，2019).随着周期的增加，鄂尔多斯盆地的四周逐渐连成一片呈低速异常的分布特征，内部整体上呈现明显的高速异常，这与盆地内部相对稳定、地震活动比较弱、断裂比较少，而周围分布一系列断陷盆地的特征相吻合.

周期为90～150 s的图像对深度在82～330 km范围的S波速度比较敏感，主要反映了深部岩石圈的结构特征.鄂尔多斯盆地下方的高速异常能够延续至200 km，说明盆地下方具有较厚的岩石圈结构，这和已有的研究成果吻合较好(An and Shi，2006；Chen et al.，2009；Chen，2010；李多等，2012；Tian and Zhao，2013；于勇等，2016；陈兆辉等，2018；高翔等，2018).鄂尔多斯盆地在岩石圈范围内以北纬38°为分界线，北部是低速异常，南部是高速异常.随着周期的增加，鄂尔多斯盆地在岩石圈范围分界线两侧的速度差异变得越来越明显，这和已有研究具有较好的一致性(Tian et al.，2009；高翔等，2018).

4 讨论

相速度的成像结果表明，研究区的主要特征是鄂尔多斯盆地逐渐的被低速异常包围，整个华北克拉通中部为显著的低速异常.其次以北纬38°为分界线，山西断陷带的南北部速度差异在下地壳和上地幔比较明显，而鄂尔多斯盆地的南北部速度差异在深度200 km附近比较明显.

研究发现整个研究区以北纬38°为分界线的南北部从浅部沉积层到深部岩石圈范围内均存在明显的速度差异的分布特征.分界线的北部有银川盆地、河套盆地、大同火山区等一系列的断陷盆地，上地壳相速度的分布特征主要受到盆地沉积层的影响.方盛明等(2009)利用深地震反射探测剖面获得银川地堑内新生代沉积最厚约9 km.Wang等(2017)利用接收函数研究表明，分界线北部的沉积层比南部的厚，其中河套盆地沉积层的厚度超过5km.对于山西断陷带在下地壳和上地幔范围内南北部差异的成因，周铭等(2016)推测山西断陷带两侧的构造应力和构造运动方式不同导致地幔物质上涌.鄂尔多斯盆地在岩石圈范围内的南北部差异，Liu等(2004)认为鄂尔多斯盆地北部上地幔存在的低速体是由于青藏高原上地幔物质运移造成的，已有研究(Tian et al.，2011；Dong et al.，2014)认为是太平板块俯冲引起地幔热物质上涌造成的.华北克拉通下方的岩石圈整体上从西到东逐渐减薄(Chen，2010；朱日祥等，2011；黄方等，2015；于勇等，2016)，这说明华北克拉通的岩石圈正在遭受不同程度的破坏.对于鄂尔多斯地区的构造应力场，目前认为鄂尔多斯块体主要受到青藏高原的挤压及太平洋板块俯冲的共同作用.前人研究(Huang et al.，2008；常利军等，2011；于勇等，2016；Guo and Chen，2017)在鄂尔多斯盆地南部存在一个狭隘的通道，来自青藏高原的软流圈通过绕流进入鄂尔多斯盆地的东南部，对鄂尔多斯盆地东南部的岩石圈进行侵蚀和改造.目前我们对这一说法持有怀疑的态度，假设来自青藏高原软流圈绕流观点成立，软流圈经过鄂尔多斯盆地南部时应该对其破坏更加严重，然而成像结果表明鄂尔多斯盆地南部的岩石圈比北部的厚.已有研究表明大同火山群具有大范围的新生代幔源岩浆岩出露(Liu et al.，2004；Xu et al.，2005)，大同火山区从下地壳开始不断的向周围延伸，随后整个华北克拉通中部为显著的低速异常，相比于来自青藏高原的软流圈绕流到鄂尔多斯盆地东南部的观点，我们更倾向于太平洋板块俯冲引起深部热物质上涌和运移这一观点.Lei(2012)成像结果表明大同火山和渤海海底的低速异常在200 km以下相连呈一个Y形结构且延伸至下地幔，认为大同火山尽管与太平洋板块深俯冲形成的“大地幔楔结构与动力学”(Lei and Zhao，2005，2006；Huang and Zhao，2006；Lei et al.，2013；Zhang et al.，2020)有关，但更可能的是起源于下地幔的地幔柱(Lei，2012).本文结果也揭示出大同火山下方从地壳到上地幔存在一个显著的柱状低速异常，因此，我们认同大同火山更有可能起源于下地幔的地幔柱这一观点.

5 结论

本文使用基于程函方程的面波层析成像的方法，利用流动台站和固定台记录到的地震面波数据，得到了鄂尔多斯及邻区瑞利波相速度的分布图.相速度分布图很清晰的反映出研究区范围内的地壳和上地幔的速度结构的精细特征，与地质特征吻合较好.

短周期时，瑞利波相速度的高低速异常的分布特征与地形地貌很好的对应.中等周期时，鄂尔多斯盆地的四周逐渐的被低速异常包围，这和鄂尔多斯盆地周围分布一系列的断陷盆地相吻合.长周期时，鄂尔多斯盆地和山西断陷带以北纬38°为分界线南北部的岩石圈厚度差异比较明显.

大同火山区从下地壳到上地幔不断地延伸和扩展随后整个华北克拉通中部为显著的低速异常，其次整个研究区以北纬38°为界从浅部沉积层到深部岩石圈的范围内南北部的结构特征存在明显的差异，我们认为主要由于太平洋板块俯冲的作用引起热物质上涌和运移从而导致鄂尔多斯及邻区的岩石圈遭到不同程度的破坏.

致谢感谢中国地震局地球物理研究所“中国地震科学探测台阵数据中心”为本文提供的地震数据资料；感谢两名审稿专家的评阅.