地震透射层析成像

1. 地震透射层析成像

1.地震透射层析成像的基本原理
地震波透射层析成像是利用地震波对于地质体的透射投影，来重新构成地质体内部地震波速度的分布形态。根据地震波速度与地质体的对应关系，进行岩土的分类和评价。
地震波速度和岩体的特性，一般都具有较好的对应关系。致密完整的岩体地震波速度较高，而疏松破碎的岩体地震波速度偏低。当地下有低速介质存在时(视为异常体)，地震波穿透这些低速介质将产生时间差。不能根据一条射线的时间差来判断低速介质的具体位置，因为它的位置可能在射线的任何一处。如果有多条射线在同一介质中穿透，并且这多条射线构成相互交叉的致密射线穿透网络，就能对低速介质在空间位置上有较强的限定。层析成像就是应用适当的反演计算方法构制速度图像，从而获得低速介质的分布位置。
2.反演计算与图像生成
利用射线理论，首先把被测区域分成若干个面积相等的成像单元，实现透视空间的离散化。设被测区域共有m个成像单元，Uj为第j个成像单元的慢度，对于每条地震波射线可写成:

物探数字信号分析与处理技术

式中i=1，2，…，n;j=1，2，…，m;n为射线总条数;m为成像单元个数;ti为第i条射线的旅行时间。
将每个成像单元的Uj(x，y)视为常数，则式(13-2-1)可写成级数形式:

物探数字信号分析与处理技术

lij为第i条射线通过第j个成像单元内的长度。对于每条射线都可以写成一条射线方程。有n条射线，便得到n个射线方程，因此式(13-2-2)实际是一个线性方程组:

物探数字信号分析与处理技术

求解式(13-2-3)可得到每个成像单元内的地震波慢度值Uj，其倒数即为地震波速度值vj，然后再采用平滑插值数字技术，绘制地震波速度等值线图。
对方程组(13-2-3)的求解方法有多种，如联合迭代法(SIRT)、代数重建法(ART)、共轭梯度法等。联合迭代法在解大型稀疏线性方程组方面具有收敛速度快及占用内存少等优点。
代数重建法的基本思想是先给出每个单元之间的初始值Uj，然后将所得到的投影值残差逐一沿其射线方向均匀地反投影回去，不断地校正Uj值，直到满足要求为止。若用数学表示上述过程，则有 初始值

物探数字信号分析与处理技术


图13-2-1 区域划分及射线分布示意

式中U(k+1)j为第j个单元第k+1次迭代的慢度值;U(k)j为第j个单元内第k次迭代的慢度值，λ为阻尼因子(0<λ<2)。
每次迭代的新迭代值U(k)j要与前一次迭代结果U(k-1)j进行比较，如果差值小于给定的精度要求E，即

物探数字信号分析与处理技术

说明已经收敛，否则还要做下一轮迭代，直到满足(13-2-5)式为止。
由于给出的收敛标准在实际工作的反演迭代过程中有时很难达到，为提高效率，可在反演时给出最大迭代次数kmax，用以双重制约程序的运行。
实际运行反演程序时，首先调入激发点、接收点坐标及地震波旅行时间的数据文件，并根据菜单提示，先后输入被测区域x和y方向的总长度及x和y方向的单元长度，并输入单元内地震波的最大速度值;其次输入收敛精度E、阻尼系数λ、最大迭代次数kmax。输入完成后进行反演计算，反演迭代完成后，可得到每个单元的速度值。启动绘制等值线程序，并输入相应控制参数，即可在Autocad的支持下输出被测区域内的地震波速度等值线图。
3.应用实例
(1)对岩体进行工程地质评价
黄河大柳树水利枢纽坝址区经历多次构造运动，断层裂隙极为发育，岩体被切割成大小不等的块体，完整性较差，加之该区处于干旱的沙漠气候，风化作用强烈。为了对该区岩体进行工程地质评价，利用现有的平硐进行了穿透测试。测试时采用相邻平硐之间相互激发与接收，激发点距为5m，接收点距为2m，使用炸药震源。利用ES-1210F工程地震仪接收纵波旅行时，这样在被测区域形成了纵横交叉的穿透射线。反演计算时，将被测区划分成2m×2m的规则成像单元。利用射线在成像单元内所经过的路径和相应的旅行时，建立线性方程组。反演采用代数重建法，得到每个成像单元的地震波速度值，而后平滑插值，并绘制成地震波速度等值线图。图13-2-2是PD304～PD306的测试成果。

图13-2-2 PD304～PD306层析成像波速图

从图可见，边坡山体内部随着风化程度的减弱地震波速度增加，另外，岩体的局部构造和结构对波速的影响也很突出，在断层破碎带、裂隙密集处。地震波速度明显下降，如PD304洞底处，在波速图像中都有明显反映。
(2)在工程检测中的应用

图13-2-3 测试点布置示意


图13-2-4 成像结果等值线图


图13-2-5 成像结果立体图

首都机场高速公路桥某桥台，台身长34.5m，宽4m。在施工中因中途更换灰浆，怀疑局部可能存在缺陷。为查明是否确实存在隐患，采用声波对穿方法，在桥台两侧12m段内布设了测试系统。接收采用高灵敏加速度计，震源采用专用手锤，激发和接收点距为1m。在指定的12m段内，共计121组数据，使得射线基本覆盖了被测区域，见图13-2-3。反演采用最优化准则和联合迭代重建技术(SIRT)。计算结果分别用等值线图和立体图示于图13-2-4和图13-2-5。
显而易见，图中有一个明显低速区，范围不大。但由于其应力波速度值仍较高(≥3000m/s)，根据《钢筋混凝土设计规范》中关于混凝土强度与弹性模量之间的关系得知，该低速区的混凝土强度≥200号，满足设计要求，工程质量合格。经调查，该低速区是由于施工存在振捣器漏捣现象。在清理台身周围的堆积残留物和模板时，发现在低速区的地方，混凝土表面呈蜂窝状，正是漏捣所致。

2. 地震折射层析成像

设地表排列方向为X轴，深度方向为Z轴，速度为v(x，z)，则二维情况下的几何地震学的基本方程为: 折射层析就是求解上述时间场方程，重建时间场，确定折射界面。
为求解(13-4-1)式，可近似给出速度曲线v(z)(由折射资料和钻孔资料确定)，并将微分变成差分，将xz剖面划分成网格状，如图13-4-1所示。沿X方向每隔Δx作一条平行于Z轴的直线，Δx可等于道间距，也可为道间距的1/k。
设沿X方向的编号为i，i=1，2，…n;xi+1-xi=Δx;沿深度方向每隔Δz作一条平行于x轴的直线，编号为j，j=1，2，…m;zj+1-zj=Δz，在深度为zj、水平位置为xi的结点处记为c(i，j)。
为重建时间场，需要求出每个结点上的弹性波到达时间t(i，j)，为此将(13-4-1)式写成差分格式:

物探数字信号分析与处理技术

式中v(i，j)表示结点c(i，j)处的速度值。对Δt/Δx按中心差分格式计算:

物探数字信号分析与处理技术

则由(13-4-2)式可解得

物探数字信号分析与处理技术


图13-4-1 地震折射层析剖分示意图

式(13-4-4)即一阶差分方程，可由龙格—库塔法求解。以地面(j=0)实际观测到的时距曲线t(i，0)在各结点的弹性波初至时间作为初始条件，经龙格-库塔递推公式可依次求出zj深度上各结点弹性波的到达时间t(i，j)。
根据折射波相遇时距曲线，地下每个结点可求出两个时间值t1(i，j)和t2(i，j)。
如果某结点正好落在折射界面上，由互换原理有:

物探数字信号分析与处理技术

TR为互换时间，当结点C(i，j)不在折射界面上时，(13-4-5)式不成立，则

物探数字信号分析与处理技术

若结点上位于折射界面之上，则上式中的差值时间e(i，j)0。设在结点C(i，j)处，e(i，j)>0在结点C(i，j-1)处，e(i，j)<0，则xi处的界面必在深度zj和zj-1之间通过。xi处的界面深度可用内插公式求出

物探数字信号分析与处理技术

对于不同的xi点，重复上述计算，就可求出整个剖面的界面深度。

3. 天然地震层析成像

王小凤等（1997）用赵大鹏等（1991）开发的一种新的地震波速图像处理程序，对辽宁省地震台网观测的数据进行处理。得出海城震区地壳在深度分别为2、9、17、26km和35km的P波速度图像。显示下辽河盆地及其邻区的地壳三维速度结构具有如下特征：
（1）在平面上天然地震震中主要沿郯庐断裂带呈北东向展布，从三维空间上看，它们主要沿郯庐断裂带的界面和地壳12～25km的低速层展布，于12～22km更为集中。由于反映郯庐断裂带地壳深部产状近直立且微向东倾，直切地壳31km深。该部位在浅部相当于下辽河盆地东部凹陷与东部凸起分界带。
（2）沿郯庐断裂营口-海城一带，在地壳深部有P波低速体集中分布带，但具体展布方式随深度有变化。在地壳2km深处，于下辽河盆地和渤海湾沉降区呈现速度为3.70km／s左右的低速区，反映了山地和盆地地壳介质的差异。低速体局部轴向近SN向。在9km深P波速度分布与2km处大体一致，于辽东山地和下辽河盆地过渡地带的营口、盖县还出现了一个仅南北向分布，P波速度为5.80km／s左右的低速区。海城震中及其附近地壳17km深处的介质P波速度相对其周围偏低，这个低速区呈北西向分布，该低速区中心相对于9km深处低速区中心向东偏移，范围也有增大。特别是低速区中心出现了P波速度为6.10～6.15km／s，比上层介质速度偏低的椭圆形低速区，也呈北西向展布，其长轴大约为50km，短轴大约为20km。海城震中及其附近地壳26km深处的介质P波速度相对其周围介质速度低，但比上层介质速度高，其速度为6.35～6.57km／s，这个低速区呈近东西向分布，低速区中心比17km深处低速区中心进一步向东移动，范围更加扩大。从总体上看，P波低速体主要沿郯庐断裂带东界隆起带边缘分布，35km处表现十分典型。低速体显示东西向及南北间轴向复合形态，但以东西向延伸为主。向深部低速体出现的范围有加大的趋势。P波低速体主要沿郯庐断裂带成带展布，越向深部显示越加明显，东界断裂带低速体分布相对更加集中，为现今强烈活动断裂。
（3）从海城—营口地壳三维速度结构变化来看，随着地壳深度从9km至35km，低速区范围不断增大，逐步向东扩展。
（4）从内蒙古乌珠穆沁旗-辽宁东沟地学断面资料分析，以上由天然地震层析成像所给出的P波低速带位置与地壳拉伸减薄，莫霍面与软流层顶界面上隆的部位相吻合，并伴随上地幔铁镁物质上涌，因此郯庐断裂带可为上地幔上涌的通道。

4. 地震层析成像的介绍

英文：seismic tomography释文：仿效医学上用x射线对人体内部组织结构进行逐层剖析成像的原理，利用地震波在不同方向投射的波场信息，对地下介质内部精细结构（例如速度、衰减系数、反射系数等的分布）进行成像。方法与反演类似，一般也将地下划分为网格，对每个单元给定要成像的参数值，计算走时或振幅，比较计算结果与实测数据，迭代进行至两者之间的误差达到要求为止。具体算法也有多种，例如代数重构、截断QR分解法、正交分解最小二乘法、截断奇异值分解法、散射层析成像法等。地震层析成像也有反射、绕射和透射层析成像之分，后者要采取井间或地面一井间观测。

5. 地震层析成像的定义

地震层析成像（seismictomography）是指利用大量地震观测数据反演研究区域三维结构的一种方法。其原理类似于医学上的CT，但地震层析成像比医学上的CT技术更复杂。大量数据以及其他许多不定因素，包括存在多种数据误差、解的不唯一性在内的地球内部成像问题。

6. 地震层析成像反演方法

地震观测数据与地球参数之间存在着积分关系。例如，地震波速度和地球的衰减性质就与实测的走时及地震波振幅有关，这种关系可以由沿射线路径的线积分获得。同样，高频电磁波的传输损耗与地球的电学性质有关系。因此，地球物理学的层析技术就是对这些积分关系进行反演，以获得对射线穿经的某些空间区域内波速场v（x，y）或电磁耗损系数场α（x，y）的估计（图4-18）。

图4-18 地球物理层析反演技术方法框图

本节介绍矩阵反演与傅氏变换方法原理，其他的地震层析成像反演方法可参考有关文献。
（一）矩阵反演
广义矩阵反演、阻尼最小二乘或线性规划等多种方法都容易应用于这一问题。把待研究的区域划分成许多单元的网格阵，并假定在每一个单元内v（x，y）或α（x，y）均为常数。以地震走时方程为例，其线积分可近似表达为

固体地球物理学概论

式中：∆sj是射线穿经第j个单元的距离；vj是第j个单元内的地震波速，求和是对第k条射线实际穿过单元来实施。在具体应用时，走时方程通常是依照求解初模型的速度扰动量来建立，即

固体地球物理学概论

这里

固体地球物理学概论

即为第j个单元内慢度的扰动量。这样处理是因为特定资料的数量还不足以唯一的确定单元网阵内所有的参数值。换言之，即未知量的数目多于线性独立的方程数。在这种情况下对方程组求解时还必须给以附加的约束，而且通常的求解过程实际上是搜索这样一个解，即它不仅要满足数据资料，还应“靠近”预置的初始模型（如对所有j单元的  值）。
矩阵方程的形式为
Y=AX
式中：Y为走时扰动值∆tk（k=1，2，…，m）组成的m维列向量；X为慢度扰动值∆Pj（j=1，2，…，n）组成的n维列向量；A是∆s（kmax×jmax）的矩阵，kmax为穿经研究区域的射线总数，jmax为单元总数。A为一种比较稀疏的矩阵（含有大量的零元素），因为在整个研究区域内，任何一条射线都只能穿经一小部分单元。
要构成A矩阵就必须知道射线从震源到接收点的路径，而这又取决于波速场，所以这一问题的答案恰是待求的解。在实际反演时，射线路径将通过初始（设想的）模型进行追踪，并由此构制A矩阵。方程Y=AX可由一些矩阵反演程序解出。X中的速度扰动量，在反演计算时通过适当的阻尼处理使其保持为小量，由此引起射线路径的偏差也假定为很小量，且忽略不计，一旦得到新的修订速度模型，地震射线就可通过该模型进行追踪并构制新的A矩阵。这一过程可能需要重复多次才能大大降低观测与计算走时的残差。
对于一个成功的层析反演来说，地震射线的角度覆盖范围必须尽可能的宽，而在大多数的地球物理应用中，震源-接收点的位置分布均受到严重的限制。这就不可避免地造成了所得图像的混乱和含糊。图4-19给出了一个相当极端的例子。通过图4-19（a）所示的简单速度模型，由10个爆炸点的位置到14个接收点的位置，对140条射线进行了追踪。
研究区宽为8个单元，长为12个单元，A是一个140×96的矩阵，图4-19（b）是通过广义逆矩阵反演方法得到的映像。图4-19（c）是只利用一个爆炸点对走时进行反演得到的图像。对于所得图像在射线路径方向上出现的轮廓不清晰是本例中角度覆盖范围极其狭窄的必然结果。

图4-19 用广义逆矩阵反演方法获得的速度场层析摄影图

（a）该速度场中有一方形异常区，其余部分为均质区，两者的速度差为10%；（b）反演结果；（c）仅用一个爆炸点发出14条射线所得反演结果
将矩阵反演方法用于层析成像的最大缺点是计算工作量太大。如果研究区域内的单元数为n，那么对矩阵求逆的运算是n3量级。所有10×10单元将需要106时间单位，而100×100单元则需要1012时间单位。若每一步运算用1μs，那么100×100单元求逆就要花277h，而200×200单元则需花掉两年的时光。在许多实际问题中，100×100单元的模型并非不现实，因此，寻求其他运算速度更快的计算方法是至关重要的。
（二）傅氏变换法
傅氏变换法是基于投影切片定理发展起来的，该定理表明：“在角度为θ时一个投影的一维傅氏变换，即为其原始客体的二维傅氏变换沿同一角度的切面值（Mersereau et al.，1974）。”图4-20（a）是一个速度场的分布图。在一均匀速度场中有一正方形的高速度区（由于计算上的原因，这个正方形的四边在图中略做了平滑，但这并不妨碍观察其主要特征）。这个速度场的投影就是许多平行射线以一固定角度穿经该场的走时。当射线平行于X轴或者Y轴时，其走时异常呈矩形，而当射线对X或Y轴以45°角入射时，其投影的走时异常便呈三角形。图4-20（b）是图4-20（a）的二维傅氏变换。根据投影切面定理，穿经图4-20（b）的平面即为图4-20（a）以同一角度投影的一维傅氏变换。熟悉一维傅氏变换的人不难发现，沿图4-20（b）中A—A和B—B线的切面将分别是矩形和三角形的傅氏变换。

图4-20 投影切片定理图示

（b）是（a）的二维傅氏变换，穿越（b）的平面。
例如A—A或B—B，是穿经（a）以同一角度投影的一维傅氏变换

图4-21 傅氏变换平行线束沿X方向的投影示意图

设平行的射线束沿X1方向，则其投影值为（图4-21）

固体地球物理学概论

式中：P表示走时或振幅；f（X1，X2）为待反演的介质特性。PX1（X2）的傅氏变换为

固体地球物理学概论

但f（X1，X2）的二维傅氏变换为

固体地球物理学概论

将式（4-49）代入式（4-50）并与式（4-51）比较得
PX1（k2）=F（k1，k2）k1
以上证明了切片投影定理，事实上，它可以推广到任意角度θ的投影上，如图4-22所示。

图4-22 切片投影定理推广到任意角度的投影示意图

上述讨论说明，原则上可以建立一个快速的反投影法。其基本思路如下：对于各个不同角度的投影函数的一维傅氏变换，可以构造出二维傅氏变换函数，再对此函数进行二维反傅氏变换，就可求得速度场分布。
这种方法的优点是运算速度快，而且还可以利用傅氏谱的对称特点减半数据量。该方法的缺点是在频率域内做内插时，高频损失较严重，这会影响到成像的分辨率，并且，内插计算也比较费时。针对地震学非完全投影情况下，可采用有限角投影重建的办法或采用扇形滤波器避开观测点缺失区的办法（Menke，1985），这却又会降低图像的分辨率，但不失基本形态。另外对射线的直线性及震源—接收器系统结构形态要求要严格限制。
其他地震层析成像的反演方法主要有褶积方法或滤波反投影法（Andersen et al.，1984；刘福田等，1989）、代数重建法（Gordon et al.，1970；刘福田，1988）以及波形与走时联合反演（Dziewonski et al.，1993）等，不做进一步讲述。

7. 地震层析成像的理论基础

（一）地震层析成像的数学理论
1917年，奥地利数学家拉冬（J.Radon）对于由投影重建图像的思想首次做出了严格的数学表达，解决了分布函数与投影之间的变换关系问题。
1.拉冬变换
在二维域中存在一未知的连续函数f（x，y），令Ox[y]坐标系逆时针旋转θ角，形成O[vu]坐标系。将f（x，y）沿平行于u轴方向的射线Li做线积分，并设其为Pθ（ti），经变换后得Pθ（ti）为f（x，y）在角度为θ时沿射线Li的投影值，沿众射线Li（i=1，2，…）投影，构成投影函数

固体地球物理学概论

式（4-43）即为拉冬正变换，改变θ角得到一系列投影函数Pθ（t）。
将投影函数与某一核函数褶积求得褶积后的投影函数Gθ（t），再做反投影求得目标函数f（x，y），即拉冬逆变换。投影成像应采用滤波反投影形成的拉冬逆变换为

固体地球物理学概论

其中s=ycosθ—xsinθ。实际应用时，要将式（4-44）离散化、有限化。
2.傅里叶（Fourier）投影定量
设目标函数为f（x，y），投影函数为P（r，θ），并满足拉冬变换，求得投影函数的一维傅里叶变换P*（u，θ）和二维傅里叶变换f*（u，θ），而P*（u，θ）=f*（u，θ）式为傅里叶投影定理，据此求取研究区的速度场。
（二）地震层析成像的物理基础——射线理论
射线方程和界面条件（Cerveny et al.，1977；刘福田等，1989）：假定地球是各向同性、完全弹性的成层介质，在利用走时重建速度图像时，取地震波的高频近似解，且将地震震源视作点源。设波速为v，则由震源i至接收点j的走时可写成

固体地球物理学概论

式中：Lij为射线路径。在球面坐标系中，路径的微分方程为

固体地球物理学概论

式中P=（Pr，Pθ，Pφ）为慢度向量：

固体地球物理学概论

式中：cosγr、cosγθ、cosγφ为射线的方向余弦，有cos2γr+cos2γθ+cos2γφ=1。
式（4-44）是在波速为连续函数假定下导出的，地球内部存在界面，则波入射后产生反射和折射。设两界面之间的地层内部是连续的，只要能够确定界面条件，就可在分层介质中进行射线追踪。界面条件可用相位匹配法导出，在球坐标系中的折射波有

固体地球物理学概论

式中：角标i表示入射波，j表示折射波；v1和v2分别为入射介质和折射介质的波速；sgn（）为符号函数。
射线宽度（Cerveny et al.，1983）：根据式（4-45）～式（4-48），不难得到地震波的走时。这就意味着，已经假定式（4-45）是在三维空间上的射线上求积分；理论上只有当无限高频时才能认为存在无限窄的射线（波长λ→0），而射线解对介质是“绝对”分辨的。由于地震波频率总是有限的，因此分辨能力也是有限的。为解决这个矛盾，在地震波的高频近似解中需要对标准射线方法进行修正。例如，可以把射线视作中心射线为费马射线的高斯射线，其宽度可用距中心射线的距离d（s）来度量。显然，射线宽度d（s）将直接影响地震层析成像的分辨率。根据惠更斯原理，一个波前的每个面元可看作一个产生球面波的次级扰动中心，且以后任意时刻波前的位置是所有这些子波的包络面。按照惠更斯菲涅尔原理，只当OAP—ODP≤λ/4时，次生波在P波处才产生相长干涉，以致A处的结构要影响到P处观测到的波场。据此，得到路径长度为L射线的最大宽度dm（s）近似为

固体地球物理学概论

例如，当λ=10km时，对长度为3000km和10000km的射线，其最大宽度dm（s）各约为61km和112km；当λ=5km时，对长度为100km和500km的射线，其最大宽度各约为8km和18km。这表明，远震和区域地震数据其水平分辨率是＋分不同的。这在地震层析成像中须予以重视。

8. 体波地震层析成像

一、层析成像反演方法
假定地球是各向同性的、完全弹性的分层介质,地震波视作点源,并取地震波的高频近似解。设波速为v,则由震源i至接受点j的走时可写成

新疆阿尔泰—天山地学断面地质地球物理综合探测和研究

式中Lij为射线路径,由路径的微分方程给出(Cerveny, V.et al, 1977)。用三维空间中的非均匀网各点的速度值描述介质的速度函数,并且允许存在介质参数的间断面(朱天飞等,1982;刘福田等,1984)。
由9.1式可知,由于走时取决于沿射线路径的积分,而路径同介质速度有关,因此,利用9.1式反演速度将使成像问题高度非线性。 我们可按下述方式线性化

新疆阿尔泰—天山地学断面地质地球物理综合探测和研究

式中 为观测走时; 是参考速度为v0时的理论走时, 分别为第i个地震的深度、纬度、经度和发震时间。
9.2式的第一项为速度相对于v0有小扰动δv时所引起的走时变化,第二项代表震源时空参数变化对走时的影响。据此,利用刘福田等(1989)提出的方法进行层析成像。
二、资料与初始模型
1.资料
用于本书方法的P波到时资料有两种,一种是研究区域以外的远震到时资料,另一种是研究区域以内区域到时资料。 由于远震资料的入射角范围相当窄(一般只在±30°范围内),它只能用于分辨尺度相当大的结构。 它的水平分辨可以用增加台站个数,减小台站间距来改善,但对地壳上部的垂直分辨仍是不够理想的。 与远震资料相比,区域地震到时资料的优点在于入射角范围可以在很大的范围内变化,因此对地壳上部的三维结构可以取得良好的分辨。 本项研究重点关注的是新疆北部地区地壳速度结构,因此采用区域地震到时资料。 一般来说,地震射线(地震波)穿透地下介质的深度与震中到接收台站的距离成正比,而且研究区域内的地震震源深度多在10 ～35km之间,如果研究区域选择过小,通过下地壳的地震波将很少,我们将无法能获得下地壳地震波速度分布图像。 为此,我们将研究区域扩大为北纬39°～50°和东经77°～91.5°,它包括中国新疆和吉尔吉斯斯坦等。 研究区域内有55个地震台站(图9-1),经过计算残差,最后选用了1984～1997年间发生在研究区域内的1049个地震(图9-2),每个地震一般不少于10个地震台站记录。 共用了13094个P波到时数据,到时误差小于0.3s。

图9-1 研究区域内的地震台站位置及网格划分示意图


图9-2 用于层析成像的地震震中分布图

2.初始模型
初始模型的选择按照以下原则:首先,根据地震台站和震中的分布,以及我们所关心的区域进行划分,保证所划分的网格有足够的射线通过,网格的尺寸是不均匀的,水平方向的最大尺寸为3°×2°,最小尺寸为0.7°×0.7°(图9-2),垂直方向的最大尺寸为60km,最小尺寸3km(表9-1)。 第二,不同深度处的平均速度值的选取是参考新疆地震台网观测报告中提供的速度模型,100km深度以下的速度参考Herrin(1968)模型。 第三,深度的划分主要考虑不同层内有足够的射线数和研究的需要,莫霍面的平均深度则根据用 与Pn实际观测资料所做的时-距曲线而确定。 研究区域内总共有2340个节点,即有2340个待求系数。 各个深度处的参考速度值由表9-1给出,并假定从一个深度到下一个深度,速度按线性变化。
表9-1 新疆北部地区的初始速度值


三、三维速度结构重建
图版9.Ⅰ是本书得到的不同深度处的速度分布图像。 其速度分布如色标所示,参考速度值为零,红色到绿色表示低速区,蓝色到紫色表示高速区。 由于受到地震台站和地震分布的限制,研究区域内的左上角、右下角等部分地区没有或很少有射线通过,所以未对这部分地区进行推论。
1.地壳上地幔速度图像
图版9.Ⅰ(a)是3km深度的速度图像,参考速度为5.8km/s。 该深度的天山与准噶尔盆地相衔接的山前地带和天山与塔里木盆地相衔接的山前地带为低速区,速度值在4.6～5.6km/s之间;在天山、阿尔泰山、阿拉泰山、塔尔巴哈台山等以高速为主,速度值在5.65 ～6.4km/s之间。
图版9.Ⅰ(b)是12km深度的速度图像,参考速度6.0km/s。该深度基本保持3km深度的速度分布特征。 在天山与准噶尔盆地相衔接处的山前地带仍为低速区,速度值在4.8～6.0km/s之间,但范围缩小;在天山与塔里木盆地相衔接的山前地带仍为低速区,速度值在4.8～6.0km/s之间,但范围已显缩小到新和(XIH)周围地区;在天山、阿尔泰山等山脉仍以高速异常为主,速度值在6.0～6.6km/s之间。
图版9.Ⅰ(c)是21 km深度的速度图像,参考速度6.2km/s。该深度阿尔泰山以高速为主,速度值在6.2 ～6.8km/s之间;天山为高、低速相间分布,高速占优势;在准噶尔盆地以低速为主,速度值在5.6～6.2km/s之间;塔里木河以北的塔里木盆地边缘地区为高速区,速度值在6.2 ～6.8km/s之间,它与天山之间有一条低速带,速度值在6.1km/s左右。左下角的塔里木盆地为低速区,速度值在6.0km/s左右。
图版9.Ⅰ(d)是36km深度的速度图像,参考速度6.5 km/s。 该深度的速度分布,在天山、阿尔泰山以高速为主,速度值在6.5～7.2km/s之间,其中有低速分布,速度值在6.2km/s左右;准噶尔盆地中心为低速区,速度值在6.2～6.5km/s之间,周围被高速所包围,速度值在6.5～7.4km/s之间;塔里木盆地为高速,速度值在6.5～7.4km/s之间。
图版9.Ⅰ(e)是51km深度的速度图像,参考速度7.98km/s。该深度为上地幔顶部的速度分布,天山、阿尔泰山普遍为低速,速度值在6.8～7.98km/s之间;准噶尔盆地以高速为主,速度值在7.98～8.4km/s之间,其中低速分布;塔里木盆地分为3块,中部为高速,速度值在7.98～8.4km/s之间,东经86°以东为低速,并与天山的低速融为一体,速度值在7.5km/s左右,图中左下角为低速,速度值在7.4km/s左右。
图版9.Ⅰ(f)是75 km深度的速度图像,参考速度8.05km/s。该深度在天山、塔里木盆地和准噶尔盆地均为高、低速相间速度分布。 塔里木盆地中间为低速,速度值在7.8～8.05km/s之间,东西两边以为高速为主,速度值在8.05～8.4km/s之间;准噶尔盆地中间为高速区,速度值在8.05～8.4km/s之间,南北两边为低速,速度值在7.8～8.05km/s之间;天山为小范围的高、低速相间分布。
图版9.Ⅰ(g)是100km深度的速度图像,参考速度8.15km/s。该深度在塔里木、准噶尔盆地均为高速,速度值在8.15～8.4km/s之间;天山为高低速相间分布,其中有两块大范围的低速,速度值在8.0～8.15km/s之间。
2.通过玛纳斯、尼勒克、乌苏、轮台、富蕴地震区的纵剖面速度分布图像
图版9.Ⅲ是通过玛纳斯、尼勒克、乌苏、轮台、富蕴地震区的纵剖面的速度分布图像,横坐标表示经度或纬度,纵坐标表示深度,深度范围0～220km。纵剖面的位置如图9.3(a)所示。 其速度分布如色标所示,深度划分和参考速度见表一。

图9-3 (a)纵剖面位置图,(b)玛纳斯地震区纵剖面位置图

图版9.Ⅱ中A-A’是通过轮台、玛纳斯地震区的纵剖面,深度划分和参考速度见表9-1,横坐标表示纬度。在速度图像中,北纬40.0°～42.0°对应塔里木盆地和天山南缘地区,在轮台地震区存在高速区,速度值在6.33～6.7km/s之间,在上地幔顶部存在高速区,速度值在8.15～8.33km/s之间;北纬42.0°～44.0°对应天山地区,在玛纳斯地震区存在高速区,速度值在6.33～6.5km/s之间,在上地幔中存在低速区,速度值在7.5～8.15km/s之间;北纬44.0°～46.0°对应准噶尔盆地西侧边缘,在上地幔中存在高速区,速度值在8.15～8.33km/s之间;北纬46.0°～48.0°对应准噶尔盆地北部地区,在上地幔顶部存在高速区,速度值在8.15～8.33km/s之间;北纬48.0°～50.0°对应阿尔泰山,在上地幔顶部存在低速区,速度值在7.35～7.65km/s之间。
图版9.Ⅱ中B-B’是通过轮台、玛纳斯、富蕴地震的纵剖面,深度划分和参考速度见表10.1,横坐标表示经度。 在速度图像中,东经81.2°～83.2°对应塔里木盆地北缘,在地壳中存在低速区,速度值在4.5～5.0km/s之间,在上地幔顶部存在高速区,速度值在8.25～8.5km/s之间;东经83.2°～85.5°对应天山地区,在轮台、玛纳斯地震区附近存在高速区,速度值在6.15 ～6.5km/s之间,在上地幔中存在低速区,速度值在7.65～8.0km/s之间;东经85.5°～89.0°对应准噶尔盆地,在准噶尔盆地南缘的地壳中存在低速区,速度值在5.0～5.5km/s之间;东经89.0°～90.8°对应阿尔泰山,在富蕴地震区存在高速区,速度值在6.25～6.5km/s之间,在上地幔顶部存在低速区,速度值在7.25～8.0km/s之间。
图版9.Ⅱ中C-C’是通过尼勒克、乌苏地震区的纵剖面,深度划分和参考速度见表一,横坐标表示经度。 在速度图像中,东经79.0°～84.3°对应天山地区,在东经81.0°附近的地壳中存在高速区,速度值在6.33～6.5km/s之间,在上地幔中存在低速区,速度值8.0km/s左右,在东经82.5°～84.0°的地壳中的尼勒克、乌苏地震区附近存在高速区,速度值在6.33 ～6.65km/s之间,在上地幔中存在高速区,速度值在8.15～8.35km/s之间;东经84.3°～89.0°对应准噶尔盆地,在东经85.0°附近的准噶尔盆地南缘的地壳中存在低速区,速度值在4.8～5.0km/s之间;在东经86°～88.0°的21km深度附近普遍存在低速区,速度值在6.3km/s左右;在东经88.0°～89.0°的地壳中存在高速区,速度值6.5km/s左右。
3.玛纳斯地震区及其附近地区的纵剖面速度图像
我们在北纬41°～45.2°、东经82.7°～88.3°的玛纳斯及其附近地区做了相互垂直的各十条纵剖面(图版10.Ⅲ和图版10.Ⅳ),纵剖面的位置如图10.4(b)所示,通过不同的纵剖面,可以得到深部构造随位置的变化情况。
图版9.Ⅲ是由西向东分布的近南北向的纵剖面速度图像,深度划分和参考速度见表9-1,横坐标表示纬度,纵坐标表示深度。 在不同位置的纵剖面速度图像上,北纬41.0°～42.0°由西至东对应塔里木盆地北缘和塔里木盆地,在21km深度以上,塔里木盆地北缘为低速区,速度值在5.0 ～6.0km/s之间(与山前坳陷沉积相关)。 在21km深度以下为高速区,速度值在6.0～7.0km/s之间。 北纬42.0°～44.0°对应天山地区,在北纬43.0°～44.0°的12km深度存在高速区,速度值在6.0～7.25 km /s之间。 北纬44.0°～45.0°对应准噶尔盆地的南缘,从西到东21km深度以上,普遍为低速区,速度值在4.5～6.0km/s之间(与准噶尔盆地南缘坳陷沉积相关),A5-A5’附近的速度值相对最高,速度值在5.5km/s左右。 结果表明,天山地区、塔里木和准噶尔盆地的速度结构明显不同,在天山的地壳中存在一个东西向长条状高速体,速度值在6.25～6.75km/s之间,在A5-A5’和A6-A6’附近的玛纳斯地震区速度值相对最低,速度值在6.25km/s左右。
图版9.Ⅳ是由北向南分布的近东西向的纵剖面速度分布图,深度划分和参考速度见表9-1,横坐标表示经度,纵坐标表示深度。在不同位置的纵剖面速度图像上,B1-B1’～B5 -B5’对应准噶尔盆地南缘,在12km深度以上普遍存在低速区,速度值在4.5～6.0km/s之间,低速区的厚度由B1-B1’到B-B5’由深变浅(沉积区呈箕状形态),在东经85.0°附近低速区不连续,在12km深度以下为正常的分层结构。 B6-B6’ ～B10-B10’对应天山地区,在12km深度附近有东西向的高速区存在,速度值在6.5 ～6.75km /s之间,在B6 ～B6’,B7～B7’中高速区分为东西两块,玛纳斯地震位于B7～B7’的两个高速区之间。 B8 ～B8’～B10～B10’中仅剩一块高速区。 在21km深度以下,高速结构呈现起伏变化。 结果表明,准噶尔盆地南缘和天山地区的速度结构明显不同,在天山的地壳中存在一个东西向长条状高速体,速度值在6.25～6.75km/s之间。
4.结果的可靠性分析
按照文献(刘福田等,1989),结果的可靠性分析用解的协方差矩阵和分辨矩阵来度量。协方差矩阵的对角元素的平方根即为相应解分量的标准差,本书解的标准差一般优于1.0%,最大不超过2.7%。 分辨不仅与通过网格的射线数有关,而且与射线对介质的采样方式有关。 从研究区域内地震台站分布和用于反演的地震分布,在除边角地区以外,射线对介质的采样方式十分理想。 本书用通过网格的射线数表示结果的分辨率,图版Ⅱ给出不同深度上通过网格的射线数分布图。 为了能够便于比较各层射线数的覆盖情况,给出0～100的统一标度,实际上除100km以外的各深度上,网格的射线数远大于100,甚至是它的40倍以上。 由于160km深度以下很少有射线通过,所以本书没有给出160km深度以下的速度分布图像结果。
四、讨 论
根据上述结果,我们可以在深浅构造关系、大地震与深部构造的关系、天山地震带与深部构造的关系和新疆北部大地震孕育发生的模式等方面取得许多新认识。
深浅构造关系
图版9.Ⅰ中3km和12km深度的速度分布与地表不同地质构造单元密切相关。 盆地、坳陷沉积等为低速区。 例如,从奎屯(KIT)—乌鲁木齐(WMQ)—阜康(FUK)一带为东西走向的低速区,它们与天山、准噶尔盆地相衔接的山前地带的坳陷沉积有关;在阿克苏(AKS)—新和(XIH)—库车一带为北东走向的大范围低速区,它们与天山、塔里木盆地北缘相衔接的山前地带的坳陷沉积有关,另外这一地带的地表分布有油气和沼泽(蔡仲琼,1993),这也是造成低速的原因。 而山脉以高速为主,例如天山、喀拉铁克山、阿拉泰山等。 在天山南北缘有速度梯度带存在,速度梯度带与地表大断裂带的位置与走向相对应。 例如,天山北缘的速度梯度带对应清水河断裂带。 速度梯度带两侧的不同速度分布表明断裂带两侧的构造不同,由速度梯度带可以推断,断裂带的深度至少达到12km。
图版9.Ⅰ中21km深度的速度分布与12km深度的速度分布相比有很大的不同。 天山为高低速相间的速度分布,这表明了天山的速度结构很不均匀,其中的低速区与部分熔融相关。 准噶尔盆地南缘的低速区与山前的凹陷沉积有关,这表明沉积厚度达到21km深度。 塔里木盆地北缘的低速区已不存在,这表明山前凹陷沉积深度不超过21km深度。塔里木盆地所能重建的北部阿克苏—新河—库尔勒地区为完整的高速区,这表明塔里木盆地结构的相对完整性。 在天山南北缘与塔里木盆地北缘和准噶尔盆地南缘相衔接处的速度梯度带,与地表的天山南缘断裂带的位置与走向基本相对应,但准噶尔盆地南缘的速度梯度带显得不很光滑,在准噶尔盆地的西北缘与阿尔泰山相衔接处的速度梯度带与地表的额尔齐斯断裂带和可可托海-二台断裂带的位置相对应,但走向不完全一致,这表明速度梯度带与断裂带相关,断裂带在地壳中表现为速度梯度带,据此推断,地表断裂带的深度超过21km。
图版9.Ⅰ中36km深度的速度分布与12km和21km深度的速度分布相比有很大的不同。 天山地区的速度分布表明结构很不均匀。 天山南北缘的低速区已不存在,这表明沉积的深度不超过36km。 与地表的断裂带相对应的速度梯度带也无明显的痕迹,这表明断裂带的深度不超过36km。塔里木盆地所能重建的部分的速度分布表明结构比较完整。 总体来说,该深度的速度值普遍较高,
图版9.Ⅰ中51 km深度的壳幔边界速度分布所作的推论有两种可能:一种是低速被解释成该地区出现异常壳幔边界,它可能表征存在地幔物质渗透引起的局部熔融;另一种是据此推论地壳厚度的变化(壳幔边界的起伏)。 按照后一种推论,在塔里木和准噶尔盆地,速度普遍大于7.98km/s,为上地幔顶部的速度,表明地壳厚度小于51km,是壳幔边界的隆起区。 在天山、阿尔泰山等山脉,速度普遍小于或等于7.98km/s,表明地壳厚度大于或等于51km,是壳幔边界的凹陷区或等于51km。 该深度的速度分布图像,表明了不同构造单元地壳厚度的差异,或者说,表明了不同构造单元的壳幔边界起伏情况。 与地表的天山南北缘断裂带相对应的位置,不同构造单元速度差异造成的速度梯度带不向地表断裂带那样光滑,而是参差不齐,但走向基本一致。 这说明了塔里木盆地、天山和准噶尔盆地等不同构造单元的接壤部位,在壳幔边界附近的构造是十分复杂的。 而且,在东经86°以西的天山壳幔边界范围与地表的天山范围相比明显缩短,有塔里木盆地高速体挤入天山低速体的现象。 值得注意是,在东经86°以东的天山和塔里木盆地地区有相联的大范围低速区,如果用壳幔边界的隆起或凹陷来解释,说明塔里木盆地在这一地区的壳幔边界情况与天山相同。 这样的结论与S波速度结构结果是相近的(宋仲和等,1993)。 但这与重力等其他地球物理结果不一致(殷秀华等,1998),因此,有待进一步深入研究。
图版9.Ⅰ中75km深度的速度分布似乎已看不出不同地质构造单元之间的差别,如果我们以地表的地质构造的划分来看该深度的速度图像,我们会得到:天山为高、低速相间的速度分布。塔里木盆地西侧为高速区,东侧为低速区。 准噶尔盆地以低速为主。 地表已知的各大地质构造单元的轮廓特征已不存在。
图版9.Ⅰ中100km深度的速度分布能够比较清楚的区分出不同的地质构造单元。 塔里木、准噶尔盆地为高速区。 天山为高、低速相间分布。各大地质构造单元之间的边界可以区分。该深度被认为是大陆岩石圈的平均深度(黄汲清等,1980),该深度的速度分布表明不同地质构造单元的岩石圈厚度是有差异的。
图版9.Ⅳ和图版9.Ⅴ中的低速区与准噶尔南缘的山前坳陷沉积密切相关,沉积的最大厚度不超过21km。
图版9.Ⅲ是3条通过新疆北部不同地区的纵剖面速度分布图像,它们的速度分布表明,研究区域内速度的不均匀性十分显著,这种不均匀性不仅表明了不同地质构造单元之间的速度结构差异,也表明了在同一地质构造单元的不同深度上构造的差异。 在不同的速度剖面中,强震与速度变化最剧烈的地壳(深度范围在10～30km)、壳幔边界和上地幔低速区相对应。
图版9.Ⅳ和图版9.Ⅴ中的玛纳斯及其附近地区纵剖面的速度分布图像表明,不同地质构造单元的速度结构不同,近南北向的纵剖面速度分布表明塔里木盆地、准噶尔盆地和天山的速度结构不同,天山南北缘两侧的低速区与天山山前坳陷有关,天山北缘沉积的厚度可达20km左右。 天山则普遍为高速区,地壳中有高速体存在。 近东西向的纵剖面速度分布表明,在B5-B5’纵剖面两侧的速度结构明显不同,靠准噶尔盆地一侧为低速区,与山前坳陷沉积有关,地壳中的速度以均匀分层结构为主,靠天山一侧普遍为高速区,在地壳中玛纳斯及邻近地区有高速体存在。
概括起来说,新疆北部地区地壳上地幔的速度结构总体上与地表的地质构造相关,随着深度的增加相关性减弱。 在天山南北缘存在山前坳陷沉积,厚度可达20km左右。 天山南北缘断裂带的深度超过21 km。 天山地区的构造不完整,存在近南北向构造带。 塔里木、准噶尔盆地地壳厚度小于51 km,天山地区地壳厚度大于51 km。 在天山东经86°以西的壳幔边界,受塔里木、准噶尔盆地的挤压明显缩短。 天山部分地区上地幔中存在低速异常。 也存在深部速度分布与地表地质构造不相关的地方,例如,在51 km深度,东经86°以东地区天山与塔里木盆地的速度结构相同。