地震层析成像LTI-SIRT算法研究与分析
摘要:本文研究反演二维近地表速度结构的初至波层析方法, 正演采用LTI算法,反演采用SIRT算法。经过长时间的实践过程,本文总结出LTI-SIRT算法具体实现方法与步骤，并对LTI 算法和SIRT算法进行了改进和发展。为了得到更好的成像结果,还结合平滑,插值,加权的处理手段,给出实验结果,分析实验误差,对算法作出评价。 关键词: 初至波；层析反演；LTI ；SIRT；射线追踪 1 引言
层析反演一般分两步：一是正演,即计算旅行时及射线追踪；二是反演求解，修正模型。不断地重复这两步,直到模型修正到满足一定的要求。射线追踪的算法有多种,比如最短路径法,有限差分方程法,旅行时线性插值法(LTI ,Linear Traveltime Interpolation)等.反演求解的方法有反投影法(BPT),代数重建法(ART),同时迭代重建法(SIRT),奇异值分解法(SVD),最小二乘QR分解法(LSQR)等多种算法。
各种算法都有自己的优缺点,经综合考虑,本文采用的是射线追踪精度高的LTI法以及比较容易实现的SIRT法。以下给出本文利用LTI-SIRT进行层析反演的实现方法与步骤。 2 用LTI-SIRT 算法进行层析反演的实现方法与具体步骤 第一步：初至拾取； 第二步：建立初始模型；
第三步：射线追踪,找出地震波从炮点到接收点的传播路径,并算出理论走时； 第四步：根据理论走时与实际拾取走时的差值来反演求解，对模型进行修正； 第五步：对修正后的模型进行插值处理； 第六步：对修正后的模型进行平滑；
第七步：重复上述四个步骤,直到模型修正满足一定的要求； 下面对以上步骤进行详细的说明. 2.1 初至拾取
初至拾取也就是对初至波的到达时间的记录，这是进行层析反演计算的基础数据。数据的好坏直接关系到层析成像的效果好坏。如果是进行实际资料计算时，初至拾取则为野外采集所得数据；若为理论模型的计算，则是在理论模型上进行正演计算所得的初至波旅行时。 2.2 建立初始模型
初始模型的建立，就是按照实际模型的地表起伏给定一假设的速度场。为方便起见，这个速度场可以是从地表开始以某一初始速度以相同的梯度往下递增。这一初始速度可以同先验知识给出，又或者从初至拾取中获得。炮－道比较接近时，可以把地震波看作直线传播，已知炮检距及初至时间就可以算出地表速度。初始模型的深度一定要够大，能够让射线自由地传播，不会出现遇到模型底部强迫反射的现象。 2.3 射线追踪
射线追踪的准确与否是影响层析成像的关键。鉴于日本科学家Aszkawar提出的LTI算法的高效率和高精度，本文采用该算法进行旅行时计算及射线追踪。LTI算法基于Fermat原理,即地震波沿着一条传播时间最短的路径进行传播。
该方法把模型离散成均匀的正方形单元，旅行时和射线路径的确定只与单元边界上的点有关。假设单元边界上任一点的旅行时可由该边界上相邻两个离散点的旅行时线性插值得到。如图1所示为一匀速正方形单元，单元边界平行于坐标轴，A,B是二个旅行时己知的点，要求射线穿过A,B边界到达D点的最小旅行时及射线路径。设射线从C点通过，C点旅行时可用线性内差公式由A,B二点的值表示。D点的旅行时为C点旅行时与波在C, D间直线传播时间之和。然后根据Fermat原理，就可求出C点的位置（即射线路径）和D点的最小旅行时。该方法也分成向前和向后处理，向前处理只计算各单元边界上的节点的旅行时；向后处理根据Fermat原理追踪射线路径．确定的射线路径不是单元边界上离散节点的连线，而是穿过单元边界上正好满足最小旅行时条件的那一点的线。这样，在一个单元内射线路径为直线，同时边界上的折射角随入射角连续变化。所以，该方法具有较高的精度和计算速度。
图 1 匀速正方形单元旅行时线性插值几何关系
A、B 为已知点，C是内插点，D是待求点
考虑到要反演地下速度结构，对纵向的分辨率要求要高，所以网格纵向长度要比横向长
度小。网格取得越小，分辨率越高，但计算量和数据存储量同时也会加大，所以要综合考虑计算效率以及内存因素。本文在对模型离散化时，采用矩形网格，计算节点只取网格线的交点。本文采用的网格模型如图2所示：
图 2
LTI存在的问题是，不能追踪到逆向传播的射线路径。这是因为，原LTI的计算方向只考虑了初至波正向传播的情况，而没有考虑逆向传播的情况，当模型速度变化复杂时，就不能追踪逆向传播的初至射线。此外，原LTI算得的节点旅行时并不一定都是最小的。原LTI算每列未知的节点旅行时时，算的节点先后顺序不一样，就可能算出不一样的值，而且无法确定先算哪个节点后算哪个节点。因为根据最小走时原理，应该是先算地震波先到达的节点也即旅行时小的节点，后算地震波晚到达的节点。可原LTI算法一开始无法确知地震波先到达哪个节点后到达哪个节点。为了解决这些问题，本文对原LTI作出了改进，具体步骤为：
1，向前处理：
（1） 按原LTI得到各节点的旅行时；
（2） 设置一个记录节点旅行时是否有变化的标志数组，数组元素对应每个节点，初始值
都为1；
（3） 除了炮点所在网格边上的节点，其余的节点再算一次。此时插值的边界包括把该节
点围起来的各个边界，如图2－3所示，如果边界两端点的标志都为0，或者两端点的旅行时都比该节点的大，则不再计算此边界。如果算到比原来小的，就用新算出来的值替代原来的。
图 3
（4） 比较此次计算得到的每个节点的旅行时和上一次计算得到的节点旅行时，如果不相
等，则置与该节点的标志为1，如果相等,则置为0；
（5） 重复步骤（3），（4），至到没有节点的旅行时不再变小为止。此时便能得到节点的
最小旅行时。 2， 向后处理：
在向后进行反向追踪时，也要考虑来自各个方向的射线。找出来自把该接收点围起来的各个边界射线中，到达该接收点旅行时最小的射线，确定此射线来自的边界以及边界上的插值点，把该插值点作为下一个接收点，以此类推，至到接收点为炮点，就完成了一个接收点到炮点的射线追踪。在此过程中，要记录每条射线经过每个格子的长度。 3，计算理论走时
利用每条射线经过每个格子的长度以及经过格子的慢度，便可求得每条线的理论走时。
上述的改进虽然可以解决原LTI算得的节点旅行时不是最小的问题，却降低了计算效率。在一些复杂的模型里，它的效率会是原LTI的1/5甚至更低。所以为了兼顾精度和效率的要求，可以设定循环次数。另外，在比较粗糙的初始模型初始迭代可以用原LTI来进行射线追踪，到后面反演得到的模型精度比较高之后，再用改进LTI。
此外，在实际工程当中，每炮和它的接收点对应的射线只会覆盖整个区域很小的一部分，所以为了提高效率，在对每一炮进行射线追踪时，可以只计算这些射线覆盖的区域，而不必计算整个区域。这个区域可以由炮点最左边的接收点和最右边的接收点来确定，考虑到有回波的可能，此区域可适当地扩大。
下面是计算循环次数不相同得到的射线追踪图。原模型是在坐标为（1500m ，300m）到(4000m,300m),以及从(0,700m)到（1500m,700m）有一水平的高速反射界面。在（1500m,300m）到（1500m,700m）有一坚直的相同速度的高速反射界面，水平反射界面以上是从地表初速度800m/s以每增深10m,速度增大30m/s。炮点在原点，图上所标的S的位置。每个接点相隔100m,第一个接收点在炮点右侧200m处。
