2017, Vol. 37
2. 中国地震局地球物理勘探中心，郑州市文化路75号，450002
正演试错模拟Seis83算法使用射线追踪[1]，通过反复对比观测资料和理论数据直到两者达到合适的拟合。由于需要多次重复拟合，正演模拟耗费大量时间，且不能提供模型参数不确定性、分辨率及唯一性的定量信息。
为获取地壳精细结构，Zelt等提出一种能同时获得二维速度和界面结构的射线反演方法[2-3]。该方法由模型参数化、射线追踪和阻尼最小二乘反演3部分组成。该方法在应用中有许多问题值得探讨，例如初始模型设计对结果有多大影响，反演计算如何快速收敛、提高分辨率，如何合理选取反演参数。本文通过数值试验来对以上问题进行分析和解决。
1 原理和方法 1.1 模型的参数化二维速度模型由任意变化的四边形块体组成，每层边界由任意数目和间距的边界节点以直线段形式相互连接，其节点数目和节点位置可以不同，但层边界必须从左至右穿过整个模型且不能和其他边界相交。在每一层的顶部和底部，P波速度场用任意数目和间隔的速度节点表示，顶部速度节点的数目和位置可以和底部不同，每一层速度节点的数目和位置也可以不同。在每层内，速度在水平方向和垂直方向都线性变化，层边界处的速度可以连续线性变化，也可以存在间断。若整个模型的界面速度为常数，则用单一速度点来表示，某一层的速度为常数则可以用单一的速度值来表示。
1.2 射线追踪二维射线追踪方程是一个零阶常微分方程组，可以写成以下2种形式[4]：
|
(1) |
式中，初始条件为x=x0、z=z0、θ=θ0，θ为射线和Z轴的正切角，v为波速，vx和vz分别是波速在X轴和Z轴上的偏导数。用有误差控制的龙格-库塔法求解上述方程组。当射线接近水平时，使用以x为积分变量的方程组；当射线接近垂直时，使用以z为积分变量的方程组；当射线与层边界相交时应用Snell定理。
采用变步长技术，沿射线自动改变的积分步长为：
|
(2) |
式中，α为自定义常数，步长按照局部的速度和速度梯度场调节，速度梯度大的地方步长变小；反之，步长被调大。该技术可以大大加快射线追踪的运算速度，并保证追踪的精度。
1.3 阻尼最小二乘反演在连续速度场v(x, z)中，设源点至接收点之间的射线路径为l，则走时t用积分形式表示为：
|
(3) |
在实际应用中，使用的是离散的形式：
|
(4) |
式中，li和vi分别为第i个射线段对应的路径长度和速度，所以走时是慢度的线性叠加，但其反演是一个非线性问题，需要进行线性化处理，即从初始模型出发通过泰勒展开并忽略高阶项。使用泰勒级数展开时，必须要建立初始模型，忽略高阶项后可得线性方程：
|
(5) |
式中，A为偏导数矩阵，Δm为模型改正向量，Δt为走时残差向量。偏导数矩阵包含第i个观测走时ti对第j个指定要反演参数mj的偏导数。对于每一次反演迭代，射线追踪过程中，走时残差向量和偏导数矩阵被计算出来。反演完成之后，参数修正向量被计算出来并应用到上一步的模型中。重复该过程，直到理论走时与观测走时拟合达到满意，或达到预定的终止标准为止。
从经典最小二乘解出发，为防止结果发散，在矩阵ATA的对角线元素上加常数。阻尼因子的加入说明它会沿最陡下降的方向收敛，保证其收敛性。式(5)的最小二乘解为：
|
(6) |
式中，Ct和Cm是数据和模型的先验协方差矩阵，其引入可以使解趋于稳定[5]。Cm越小越有利于迭代的收敛，但也意味着对模型的预约束就越强，或对模型的先验信息了解较多，如果盲目地减小方差可能会引入本来不存在的模型先验信息。D为全局阻尼因子，通常为了保证收敛可以先给一个较大的值，在以后迭代的每一步通过乘以一个小于1的数，以便于增大收敛速度。
模型的参数分辨率矩阵为[6]：
|
(7) |
迭代终止标准有以下2个条件：1)最终模型是在均方根走时残差RMS和无量纲指标χ2与参数分辨率两者之间取折中。一般情况下，增加参数可以减少残差，但是分辨率会降低，因此应在两者中权衡。2)尽量追踪到尽可能多的射线。如果最终模型中大部分的射线都没有被追踪到，则这个最终模型应该舍弃掉。
2 数值试验为了分析影响反演结果的因素，通过采用不同的初始模型、不同的观测走时随机扰动、不同的数据不确定性、不同的阻尼因子进行数值试验，观察反演结果与真实模型的接近程度，从而得到不同参数对反演结果和分辨率的影响。
观测系统设计如下：炮点设定为5炮，炮间距30 km，桩号从100 km到220 km，接收仪器间距2 km。“真实”模型设计为2层，如图 1所示，第一层上界面深度为0.7 km，速度为2.5 km/s；第二层上界面深度为4 km，界面上下速度相同，都为5.1 km/s；第二层底部深度为17 km，速度为6.25 km/s。利用Seis83正演追踪到的Pg波走时作为RAYINVR反演中的观测数据，如图 2所示，震相如表 1所示。
|
图 1 真实模型 Fig. 1 True model |
|
图 2 Seis83正演射线追踪 Fig. 2 The Seis83 forward ray tracing |
|
表 1 地震剖面震相 Tab. 1 Seismic phase of the profile |
在初始模型与真实模型相同的情况下，得到均方根走时残差RMS为0.034，无量纲指标χ2为0.114，直接达到迭代终止的标准。RAYINVR中射线追踪和走时拟合如图 3所示，拟合效果很好。
|
图 3 射线追踪和走时拟合 Fig. 3 The ray tracing and traveltime matching |
将第一层下速度中所有节点的速度值5.1 km/s分别修改为初速度4 km/s、2 km/s和0.1 km/s，用RAYINVR分别进行反演，得到反演结果和分辨率。阻尼因子取5，对所有反演节点的分辨率求平均值，如表 2和图 4、5和6所示。
|
图 4 初始速度改为4 km/s反演后射线追踪和走时拟合 Fig. 4 Initial velocity is changed to 4 km/s after inverse the ray tracing and traveltime matching |
|
图 5 初始速度改为2 km/s反演后射线追踪和走时拟合 Fig. 5 Initial velocity is changed to 2 km/s after inverse the ray tracing and traveltime matching |
|
图 6 初始速度改为0.1 km/s反演后射线追踪和走时拟合 Fig. 6 Initial velocity is changed to 0.1 km/s after inverse the ray tracing and traveltime matching |
|
|
表 2 不同初始速度的反演结果 Tab. 2 Inversion results with different initial velocity |
把第一层下速度中所有节点的速度值改为4 km/s进行反演，得到均方根走时残差RMS由0.404变为0.061，无量纲指标χ2由16.379变为0.379，迭代进行了3次。该行速度值可以恢复到5.01~5.20 km/s，平均分辨率为0.177。
把第一层下速度中所有节点的速度值改为2 km/s进行反演，得到均方根走时残差RMS由31.432变为0.093，无量纲指标χ2由206.115变为0.874，迭代进行了6次。反演后该行速度值可以恢复到4.91~5.41 km/s，平均分辨率为0.175。
把第一层下速度中所有节点的速度值改为0.1 km/s进行反演，得到均方根走时残差RMS由7.466变为0.113，无量纲指标χ2由5 603.453变为1.288，迭代进行了8次。反演后该行速度值可以恢复到4.68~5.26 km/s，平均分辨为0.164。
可以看出，当初始模型与真实模型差距在一定范围内时，能够基本反演到真实模型。初始模型的选择与真实模型越接近，反演结果越好。
2.2 走时随机扰动反演试验把第一层下速度值取为4 km/s，对走时数据分别加入±0.2 s、±0.4 s、±1 s的随机误差，使用相同的参数和初始模型进行反演，对所有反演的节点求分辨率的平均值，如表 3和图 7、8和9所示。
|
图 7 走时扰动为±0.2 s的射线追踪和走时拟合 Fig. 7 The ray tracing and traveltime matching with ±0.2 s traveltime random error |
|
图 8 走时扰动为±0.4 s的射线追踪和走时拟合 Fig. 8 The ray tracing and traveltime matching with ±0.4 s traveltime random error |
|
图 9 走时扰动为±1.0 s的射线追踪和走时拟合 Fig. 9 The ray tracing and traveltime matching with ±1 s traveltime random error |
|
|
表 3 不同走时扰动的反演结果 Tab. 3 Inversion results with differenttraveltime errors |
对走时数据加入±0.2 s的随机扰动，得到均方根走时残差RMS由0.426变为0.131，无量纲指标χ2由18.255变为1.716，迭代进行了4次。反演后该行速度值可以恢复到4.96~5.20 km/s，平均分辨率为0.175。
对走时数据加入±0.4 s的随机扰动，得到均方根走时残差RMS由0.461变为0.238，无量纲指标χ2由21.369变为5.690，迭代进行了4次。反演后该行速度值可以恢复到4.74~5.40 km/s，平均分辨率为0.164。
对走时数据加入±1 s的随机扰动，得到均方根走时残差RMS由0.681变为0.590，无量纲指标χ2由46.602变为34.957，迭代进行了4次。反演后该行速度值可以恢复到4.48~5.42 km/s，平均分辨率为0.114。
通过上述试验可以看出，走时随机扰动越小，反演结果越好，分辨率越高。
2.3 数据不确定性反演试验第一层下速度值取为4 km/s，数据的不确定性分别为0.2 s、0.3 s、0.5 s，使用相同的参数进行反演，对所有反演的节点求分辨率平均值，如表 4和图 10、11、12所示。
|
图 10 数据不确定性为0.2 s时的射线追踪和走时拟合 Fig. 10 The ray tracing and traveltime matching with 0.2 s identification error |
|
图 11 数据不确定性为0.3 s时的射线追踪和走时拟合 Fig. 11 The ray tracing and traveltime matching with 0.3 s identification error |
|
图 12 数据不确定性为0.5 s时的射线追踪和走时拟合 Fig. 12 The ray tracing and traveltime matching with 0.5 s identification error |
|
|
表 4 数据不确定性反演结果 Tab. 4 Inversion results with identification errors |
数据不确定性为0.2 s时，迭代进行了5次，得到均方根走时残差RMS由0.404变为0.184，无量纲指标χ2由4.095变为0.853。反演后该行速度值可以恢复到4.50~4.80 km/s，平均分辨率为0.0 543。
数据不确定性为0.3 s时，得到均方根走时残差RMS由0.404变为0.238，无量纲指标χ2由1.820变为0.988，迭代进行了4次。反演后该行速度值可以恢复到4.22~4.39 km/s，平均分辨率为0.0 268。
数据不确定性为0.5 s时，得到均方根走时残差RMS由0.404变为0.371，无量纲指标χ2由0.655变为0.404，迭代进行了1次。反演后该行速度值可以恢复到4.02~4.04 km/s，平均分辨率为0.0 104。
通过以上试验可以看出，增大数据不确定性可以降低χ2的值，但是反演结果会变差，分辨率也会降低。
2.4 阻尼因子的选取试验理论上讲，当阻尼因子逐渐减小时，反演解逐渐接近于经典最小二乘解，收敛的速度会变快；当阻尼因子逐渐增大时，收敛速度变慢，但可以保证迭代过程稳定进行。取华北克拉通超长观测剖面中三炮的数据来对阻尼因子选取进行分析。选取地震记录中能量较强、信噪比较高的初至震相Pg对研究区域进行分析，如表 5所示。
|
|
表 5 地震剖面震相 Tab. 5 Seismic phase of the profile |
试验中初始模型与真实模型相差很大，如图 13所示(方框表示深度反演点，圆点表示速度反演点)。模型分为2层、28块，60个速度参数，30个深度参数。模型深度以海平面为基准面，海平面以上为负值，以下为正值。第一层上界面深度为0.7 km的水平界面，速度范围为2.6~3.2 km/s。第二层上界面为速度间断面，深度为4 km，间断面以上的速度为5.0~5.2 km/s，以下为5.5~6.0 km/s。利用初始模型得到均方根走时残差RMS为0.346，无量纲指标χ2为1 205.34。射线追踪图中走时拟合结果很差，如图 14所示。
|
图 13 模型参数化后的初始模型 Fig. 13 The initial model after the model parameterized |
|
图 14 初始模型下射线追踪和走时拟合 Fig. 14 The ray tracing and traveltime matching under the initial model |
大的阻尼值使反演迭代次数增多，降低算法的效率，且大大降低了结果分辨率；阻尼值取得较小，可能使反演无法进行。经过对不同阻尼因子的多次尝试，得到如下结果。
1) 阻尼因子小于18时，迭代过程进行到某一步会出现计算停止的情况。这是由于选取的初始模型与真实模型差距较大，所加阻尼不足以保证迭代收敛。在理论中，应当让阻尼因子小一些，这样观测数据比阻尼对结果的影响大，但由于初始模型不准确和观测数据的影响，为保证迭代正常进行，阻尼因子不能太小。
2) 阻尼因子大于1 000时，迭代虽然可以进行，但收敛速度变慢，最后得到的χ2值会稳定到100以上。取大的阻尼值将使反演迭代加长，降低算法效率，且大大降低了结果分辨率，不满足迭代终止的要求。
3 结语RAYINVR地震走时反演算法本质上是一种线性化反演算法，从十分不同的初始模型开始反演有可能收敛到次级极值点，搜索不到全局最优解[7]。
1) 好的初始模型和观测数据能够得到好的反演结果，且能提高分辨率。初始模型的选取应当参考一维模型结果和其他先验信息。在获取观测数据时，应当尽可能准确以提高数据的质量。
2) 小的数据不确定性可以提高分辨率和改善反演结果，但在遇到观测数据较差或者χ2值过大的时候，可以适当增大识别误差。
3) 在初始迭代时阻尼因子可以给一个比较大的值，以后在迭代的每一步通过给该值乘以一个小于1的数，逐渐减小阻尼因子，向经典最小二乘靠近以便加快收敛速度。当阻尼因子过小使得迭代发散时，阻尼因子除以一个小于1的数可保证收敛性。
4) 试验数据中分辨率都很低的原因主要是射线数目较少，但不影响对比分辨率结果。
| [1] |
Ĉerveny V, Molotkov I A, Molotkov I A, et al. Ray Method in Seismology[M]. Prague: University of Karlova, 1977
(  0) |
| [2] |
Zelt C A, Smith R B. Seismic Traveltime Inversion for 2-D Crustal Velocity Structure[J]. Geophys J Int, 1992, 108(1): 16-34 DOI:10.1111/gji.1992.108.issue-1
( 0) |
| [3] |
Zelt C A, Ellis R M. Practical and Efficient Ray Tracing in Two-Dimensional Media for Rapid Traveltime and Amplitude Forward Modelling[J]. Canadian Journal of Exploration Geophysics, 1988, 24(1): 16-31
( 0) |
| [4] |
Braun J A, Smith R B. Two-Dimensional Inversion for Crustal Structure Using Refraction/Wide-Angle Reflection Data from the 1986 PASSCAL:Basin-Range Experiment[J]. EOS Trans Am Geophys Un, 1989, 70: 1206
( 0) |
| [5] |
Huang H, Spencer C, Green A. A Method for the Inversion of Refraction and Reflection Travel Times for Laterally Varying Velocity Structures[J]. Bulletin of the Seismological Society of America, 1986, 76(3): 837-846
( 0) |
| [6] |
Lutter W J, Nowack R L. Inversion for Crustal Structure Using Reflections from the PASSCAL Ouachita Experiment[J]. J Geophys Res, 1990, 95(B4): 4633-4646 DOI:10.1029/JB095iB04p04633
( 0) |
| [7] |
王夫运, 张先康, 杨卓欣, 等. 用地震走时反演长白山天池火山地区的二维地壳结构[J]. 地震学报, 2002, 24(2): 143-152 (Wang Fuyun, Zhang Xiankang, Yang Zhuoxin, et al. 2-D Crustal Structure of Changbaishan-Tianchi Volcanic Region Determined by Seismic Travel Time Inversion[J]. Earthquake Science, 2002, 24(2): 143-152)
( 0) |
2. Geophysics Exploration Center, CEA, 75 Wenhua Road, Zhengzhou 450002, China