地震层析成像doc.docVIP

地震层析成像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算每列未知的节点旅行时时，算的节点先后顺序不一样，就可能算出不一样的值，而且无法确定先算哪个节点后算哪个节点。因为根据最小走时原理，应该是先算地震波先到达的节点也即旅行时小的节点，后算地震波晚到达的