地震层析成像的正演与反演初步

地震层析成像的正演与反演初步

摘要

本文通过设立一个平行层的地球模型，初始的震源位置和发震时刻，并改变震源出射角的值，求出射线到达地面的位置，以及射线到达台站的到时，获得了正演模型得走时。并将正演结果用于反演。在反演中，本文采用了赵大鹏的反演程序1，2，反演速度结构并与设立的模型比较，得到较满意的结果。

1、引言

最初用于医学造影的成像技术自从上个世纪七八十年代引入地学后已经发展成为一项成熟的技术，越来越多地用于地球动力学，地幔对流，板块俯冲带及其演化历史，以及消亡的板块的演化历史的研究，并为板块构造理论提供有力的证据。由于到达台站的地震波的到时与地震波在所穿过的物质中的波速有关，因此，分析地震波的到时数据就可以得到地下波速结构。结合其它的地学证据，层析成像揭示出地幔由集中的上升结构与下降结构组成10。高速带通常是冷的岩石圈板块在板块的会聚边界陷入地幔的区域3，6，10，11，12；集中的低速结构通常预示着热的岩浆活动3，10，例如太平洋板块与欧亚板块碰撞形成的火山岛弧下的岩浆活动3，以及东非裂谷带下大规模的岩浆活动，导致了非洲大陆的抬升10。在对地震波的各向异性的研究中，James Wookey等8根据澳大利亚地震台站接收到的来自Tonga-Kermadec和New Hebrides俯冲带的深源地震的s波分裂，揭示出在该地区地幔中部约660km深处可能存在中部地幔分界层，阻断上下地幔的对流。随着成像解析度的提高，现在已经能反演出地球深部的速度结构和异常，追述消亡的板块的演化历史5，11，12。例如Van der Voo等10在西伯利亚1500-2800km深处发现了高速异常带，揭示了大约150-200百万年前Kular-Nera洋关闭，Mongolia-North China陆块与Omolon陆块结合的演化历史。

目前层析成像技术正向着高精确性，大数据量和适用性的方向发展，正反演数值计算方法的开发，成像方法的评价，成像结果的地学解释都是目前研究的方向。本文就是对成像的正演和反演问题进行的初步研究。

2、正演

正演是利用地质和地球物理参数定义地下构造的数值模型，计算事件在模型中的理论值。正演是提高对地下构造以及物理性质的认识的重要手段。同时在反演中也包含正演过程，因此对正演方法的研究是和有必要的。本文利用正演模型求得地震走时数据，作为反演的初始数据。

2.1模型

本文选取的计算区域是北纬35-40°，东经110-115°的矩形区域，设该地区共发生20个地震，地震发生的密度大致均匀（见图1-a），每个地震的地震参数是已知的，地震深度从240km-400km。共有25个台站均匀分散在该地区（见图1-b），台站的海拔均为0km，并假设每个台站都能清晰地记录每个地震直达波的到时，记录误差为0。由于选取的区域不大，地层在选取的区域中可近似看成是水平的，为了简便，本文采用简单的平行层模型（见图2），将400km以上分为10层，其中20km深处为Conrad面，35km深处为Moho面，35km以下分为8层，每层内速度是均匀的，其中，120km-160km层是低速层。

图1-a,震源位置 图1-b ，台站位置，充填黑色的台站为示例台站

图2，正演平行层模型 2.2方法

由于本文采用的正演模型比较简单，震源出射角度的微小变化不会导致走时的巨大差异，因此，本文采用尝试法求直达p 波的走时（见图4）。因本文选取

的区域不大，故采用 1相当于r 180π

km 的近似，记为k ，其中km r 05.6371=为地球半径。设震源的经纬度为()00,y x ，第i 个台站的位置为()i i y x ,，震中距

()()2020y y x x k

i i -+-=∆ (1)

另一方面（如图3），设震源位于第j 层，出射角为j a ，第j 层的波速为j v ，在第m 层（1≤m

p v a v a m

m

j

j ==

sin sin (2) 求得

j

j

m m v a v a s i n s i n =

(3)

2

sin 1cos ⎪⎪⎭

⎫

⎝⎛-=j j

m m v

a v a (4) 设第m 层的深度为m d ，则地震射线在第m 层的水平分量

2

sin 1sin ⎪⎪⎭

⎫

⎝

⎛-=

j j m j

j m m m v

a v v a v d l (5)

由此可求得震中距

∑

=⎪⎪⎭

⎫

⎝

⎛-=∆j

m j j

m j j

m m v

a v v a v d 1

2

'

sin 1sin

(6)

图3，snell 定律

由上式可以看出震中距只是出射角j a 的函数，由于解析解难以求出，可通过改变出射角，从0-90°，求得震中距'∆，当

ε≤∆-∆' ，001.0=ε (7) 认为该出射角j a 为真实的出射角。随后，即可求得走时

∑∑==⎪⎪⎭

⎫

⎝

⎛-==j

m j

m j j

m m

m

m

m m

v

a v v d a v d t 112

sin 1cos

(8)

图4，尝试法求射线

2.3正演模型的走时计算结果

使用上述方法即可求得模型中震源到各个台站的直达p波的走时，共获得25×20个直达p波走时数据。表1是图（1-b）中部分台站接收到的走时数据（s）。

表1，部分台站的地震到时数据

3、反演

本文的反演采用的是赵大鹏的程序，该程序被用于反演日本东北部地壳和上地幔的p波与s波的速度结构，特别是Conrad面，Moho面，以及太平洋俯冲带上的速度结构。图（5）是该程序的基本算法思路。程序的根本是步骤5，计算线性方程组，反演分辨率的高低主要在于线形方程组的解法的选择。步骤2、4也是程序的重点，因为射线路径的真实度和方程组系数的精确度直接影响到解的精度。本文将重点讨论这几个步骤。

本文计算的初始模型利用Conrad面和Moho面两个速度不连续面将所模拟的区域分成3层。程序采用两种方法来表达速度不连续面随深度的分布1，2，一是把深度作为经纬度的连续函数表达，另一种方法是使用二维节点来表达一个不连续面。当节点的深度分布已知时，不连续面上任何一点的深度可以利用周围四个节点通过线性插值的方法求得。

反演中可将地下介质分为小块体，假设每个块体中的速度均一，以块体中的一个点为代表，形成用三维节点表达的三维速度结构(见图6)。理论上，在块体