连续迭代计算方法在震波走时成像中的应用

李步青，李建楼

（1.宿州学院 信息工程学院；2.宿州学院 资源与土木工程学院，安徽 宿州 234000）

地震走时成像技术是一项系统而复杂的工程，由于它在异常体探测中具有重要的作用而不断得到发展.震波走时成像技术在震波成像中具有独特的优势，主要是因为震波走时便于测量和拾取，利用适合的路径算法能够快速、高效得到慢度分布特征，从而圈定异常体位置.但是，由于成像算法不同，反演得到的图像也不同.迭代法是解线性方程组的一个很重要的方法，适合解系数矩阵为稀疏矩阵的大型线性方程组[1,2].本文拟在Matlab软件中采用直接法和连续迭代法实现震波走时成像数值模型的反演和对比，并利用震波走时成像工程案例反映连续迭代法的成像效果，希望能为震波走时成像初学爱好者提供一定的借鉴.

1 震波走时成像基本理论

震波走时成像是典型的地球物理反演问题，涉及以下几个方面：①模型参数化，②正演（射线追踪），③反演，④解的评价[3].首先要建立模型，布置观测系统，用分块法或节点法模拟局部小异常，把物理问题转化为求解一个大型、稀疏的、且往往是病态的线性方程组问题，转化为矩阵运算问题[4]；为了得到适当的解,网格划分应该合理，以震源点间距或检波点间距划分网格大小较为合适，可提高成像精度，并非网格划分越小越好[5]；另外，要从数学上剔除病态问题,并把更多的先验信息引入反演过程[6].模型参数化之后开始正演，即根据射线路径计算理论走时.计算理论走时可以根据文[7]提出的一种方法,即假定震源到检波点路径为直线,从而可以记录每条射线穿过的单元和统计每个网格单元穿过的射线数目等，从而为反演做准备.正演是为反演做准备的，正演为反演提供了理论震波走时.反演则是利用震波的理论到时和实际到时以及对应的系数矩阵反向求解慢度值.反演方法较多，例如直接反投影法是对每条射线都按照像元中射线穿过长度占整个射线长度的比例，把走时残差分配到每个像元中，然后像元内走时残差之和除以像元内射线段总长度即为该像元慢度值，不进行迭代，不存在发散问题，但成像分辨率低[8]；代数重建法将每一条射线走时差平均地分配给通过的单元,而不考虑像元内射线的长短，是按射线依次修改像元慢度值的一类迭代算法，这种方法主要问题是收敛性不好.连续迭代法是在某一轮迭代中, 所有像元的值都用前一轮的迭代结果来修正[9,10].连续迭代法还可以考虑射线弯曲的影响,此时只需在计算理论走时值之前插入射线追踪过程即可[11].

2 震波透射走时成像数值模型及分析

图1 震波走时成像数值模型

由图1可以看出，观测范围划分为4个网格，X1、X2、X3和X4分别表示四个网格的慢度值；使用6个发射点、6个接收点，1发1收方式，合计6条射线，对应的旅行时分别为t1～t6.

每条射线的传播时间可以用下列式子表示：

写成矩阵的形式，即为A*X=T

其中，A为系数矩阵，X为慢度矩阵，T为传播时间矩阵.本例中，

假设我们不知道慢度分布情况，不妨根据经验设图1中的慢度均为100，即

令直达波理论到时为T0，则

现在给慢度一个扰动，慢度分布不再均匀，例如

按照最短走时路径追踪方法，估算得到的旅行时间为T=[180；100；160；120；198；198].

实际观测到的旅行时应和估算值差别不大，可令观测的直达波旅行时为 T=[180；110；170；130；210；200]；根据 A和T对慢度X进行直接法反演，反演结果为

其中，pinv(A)是MATLAB软件中求矩阵伪逆的函数.

直接法反演后的慢度分布与实际慢度对比及误差分析如表1：

表1 直接法反演慢度和实际慢度对比

由表1可以看出，利用直接法对慢度分布反演可以定性评价，但定量误差较大.

现在仍然使用该系数矩阵，采用连续迭代方法进行慢度反演计算.在MATLAB软件中使用迭代法对慢度逼近的过程代码如下：

100次迭代反演后的慢度分布和实际慢度对比及误差分析如表2：

表2 连续迭代反演慢度和实际慢度对比

由表1和表2对比可以看出，直接法计算最大误差为18.45%，迭代法计算最大误差为12.45%，迭代法将误差大大缩小.因此，在观测系统矩阵系数A不变的情况下，利用迭代运算方法能够有效降低误差，逼近真实慢度，提高反演精度.

3 应用案例

3.1 建立观测系统物理模型

校园内的草地中间有一条素混凝土路面，在马路的两侧按设计位置分别布置12个震源和12个检波器，形成一个震波走时观测系统.实验现场观测系统布置图2：

图2 路基观测系统布置

3.2 模型参数化

根据观测系统建立坐标系，划分慢度网格单元（即像元），模拟射线路径，根据炮点和检波点坐标计算可以得到各单元内每条射线穿过它的长度.模拟的射线路径如图3：

图3 路基探测射线模拟图

3.3 震波数据采集和处理

震波数据采集使用人工激发的锤击震源，信号采集仪器为KDZ1114-3型号矿井地质探测仪，检波器采用TZBS系列100 Hz高阻尼传感器.设置的主要采样参数为：采样频带0-4000Hz，采样间隔120μs，采样点数1024，固定增益0dB，超前采样点数0；按照激发1次、12个检波器同时接收的方式，逐点激发地震波，共可获得144道地震记录.每炮记录（12道）按放炮顺序进行拼接，拼接后的地震记录如图4：

图4 路基探测地震记录

3.4 反演结果与分析

对地震记录进行震波走时拾取和保存，然后利用连续迭代法反演，得到慢度反演值，继而得到速度反演值.速度分布图如图5：

图5 探测范围速度分布反演结果

经图5和图2对比可以看出，速度反演得到的路面位置和实际位置基本吻合，说明连续迭代方法能够应用于震波走时成像.

4 结论

根据数值模拟和实际应用结果，可以得到如下结论：

（1）利用连续迭代方法能够有效降低多元线性方程组解的误差，提高求解精度.

（2）采用连续迭代法反演基本能够对探测范围内的异常体进行定形和定位.

〔1〕 郝艳花.Jacobi迭代法与Gauss-Seidel迭代法[J].山西大同大学学报(自然科学版),2017，33（5）：3-5.

〔2〕 胡志成.Jacobi与 Gauss-Seidel迭代的比较及算法的MATLAB 实现[J].高师理科学刊,2018，38（3）：57-60.

〔3〕 和锐，杨建思，张翼.地震层析成像方法综述[J].CT理论与应用研究,2007，16（1）：35-48.

〔4〕 赵大鹏.地震层析成像及其在消减带和地震断层带成像中的应用[J].世界地震译丛,2001(2):1-8.

〔5〕 陈国金.井间层析成像质量影响因素分析[J].石油物探,1996,35（4）:18-28.

〔6〕 裴正林.井间地震层析成像的现状与进展[J].地球物理学进展,2001,16（3）:91-97.

〔7〕 段心标,金维浚.井间地震层析成像初始速度模型[J].地球物理学进展,2007,22（6）:1831-1835.

〔8〕 杨文采.地震层析成像在工程勘测中的应用[J].物探与化探,1993,17（3）:182-192.

〔9〕 刘盛东,李承华.地震走时层析成像算法与比较[J].中国矿业大学学报,2000,29（2）:211-214.

〔10〕 杨艳,秦克伟,张东，等.一种改进的近地表层析成像SIRT 算法[J].武汉大学学报(理学版),2009,55（2）:201-205.

〔11〕 雷栋,胡祥.地震层析成像方法综述[J].地震研究,2006,29（4）:418-426.