2019, Vol. 62
2. 中国石油大学(北京)油气资源与探测国家重点实验室, 北京 102249
2. China University of Petroleum(Beijing), State Key Laboratory of Petroleum Resources and Prospecting, Beijing 102249, China
黄土塬覆盖区的地震静校正问题具有两个特点.第一个特点是静校正量大.典型的黄土塬地貌特征是地表起伏、沟壑纵横.因为表层结构变化大，地表高程落差大，地震数据处理时静校正问题严重.第二个特点是静校正处理难.由于黄土塬表层覆盖巨厚黄土，近地表的低降速带厚度大，区域内潜水面较深，常规表层折射静校正方法无法达到令人满意的处理效果.黄土塬覆盖区地震静校正处理的突出难点是长波长的静校正问题(杨城增等，2016)，而层析反演静校正方法正是解决长波长静校正问题的最佳技术手段之一(张林等，2017).因此，本文探讨层析反演静校正方法解决黄土塬覆盖区静校正难题的适用性和可靠性.
在地震数据处理中，常规静校正方法包括高程静校正、模型静校正和折射静校正等(Wiggins et al., 1976; Marsden, 1993; Cox, 1999).高程静校正顾名思义仅仅校正物理测点(炮点或检波点)与固定基准面之间的高程差异，并且依据地表一致性的假设.但是在地表高程变化剧烈的情况下，地表一致性假设也许不能成立，常规高程基准面校正需由波动方程基准面校正方法来替代(杨锴等, 2007).不管什么情形，高程静校正无法适用于黄土塬覆盖区，因为该方法没有考虑低速带以及低速带内的横向速度变化状态.
对低速带横向速度变化开展调查的常规方法中包括小折射、微测井等方法，在横向上一些离散位置测定和分析近地表的地层速度，进而通过内插等算法建立近地表速度模型.根据这种低速带模型进行静校正处理的方法通常称之为模型静校正方法(查文锋和于小磊，2014).但是，小折射低速带调查方法无法适用于黄土塬覆盖区，因为该地区的地形起伏剧烈，低降速带厚度及其速度的变化较大，折射波时距曲线的拟合非常困难.即使得到了近地表速度模型，该模型也不能反映黄土塬的速度变化真实状态.微测井低速带调查方法同样不适用于黄土塬覆盖区.由于地形起伏剧烈，微测井需要加密测点，所以效率低、成本高，而且效果差.
折射静校正方法在黄土塬覆盖区也不适用.折射静校正方法利用的是地震记录上的折射初至信息，需要假设有可连续追踪的、稳定的折射界面，并且假设近地表低速带模型的横向速度变化不甚剧烈，才能有效计算静校正量.
上述三种常规静校正方法是线性静校正技术的代表，而非线性静校正技术一直也是地球物理领域的一个研究热点.Rothman(1985, 1986)利用统计学的蒙特卡罗(Monte Carlo)技术对地震道之间的互相关函数进行全局优化, 进而获取静校正量.该方法无需拾取互相关函数的峰值, 而是将互相关函数转换为概率分布, 然后对概率分布开展随机数抽样.Dahl-Jensen(1989)对最大功率谱进行全局优化, 计算基岩地区反射地震数据的静校正量.全局优化过程中使用了迭代方法, 对比上述的蒙特卡罗方法，该迭代优化方法有效地提高了计算效率.除了最大功率谱法，其他对蒙特卡罗技术的改进方法还包括模拟退火法、遗传算法等，共同的目标都是加快对全局最优解的搜寻速度.井西利等(2002)提出将模拟退火法与均匀优化设计方法相结合，林依华(2003)则提出将最大功率谱法、模拟退火法和遗传算法相结合.这些方法的综合应用期望具备三个优势，即不但全局搜索能力强，而且局部收敛速度快，同时模型的迭代更新趋于合理.Le Meur等(2011)使用模拟退火法估算强噪声数据上大量级的剩余静校正量.Abbas等(2018)提出使用模拟退火法计算复杂地表区域的剩余静校正量.
本文探讨的层析静校正方法也是一种非线性静校正技术.层析静校正(Zhu et al., 1992)采用走时层析反演方法构建近地表速度模型，然后计算静校正量.但是，与折射静校正方法不同，层析静校正方法并不强调初至波一定是折射波.由于采用了层析反演的思路，对模型空间中的任意单元进行“多次射线覆盖”，所得到的近地表速度模型具有很强的可靠性.而且层析静校正方法同步运用大量的初至信息，对每一个物理点(炮点和检波点)进行“多次走时覆盖”，最终的静校正量具有很好的统计性.正是由于上述两重“多次覆盖”特性，层析静校正方法能够解决长波长的静校正问题，是适用于黄土塬覆盖区静校正问题的有效方法之一.
本文首先对层析静校正方法的原理做出简单而又必要的回顾，给出方法实现的基本流程和技术细节.迭代反演时采用的是同步迭代重构算法(SIRT)，本文提出对同步迭代重构算法进行改进，使得层析反演的迭代过程趋于稳定.虽然层析静校正方法不受静校正量大小的制约，但是因其统计平滑的特性，当黄土塬覆盖区相邻物理点之间的高程差异较大时，仍然有剩余静校正的问题.所以，在运用层析静校正方法解决长波长静校正问题的同时，还需要运用初至波剩余静校正方法解决短波长的静校正问题.实例表明，两者结合能够有效解决黄土塬覆盖区的静校正难题.
1 层析静校正的方法原理近地表速度层析反演的基本方程是沿地震射线的走时方程，公式为
|
(1) |
式中V(x, z)为近地表介质的二维速度模型，S(x, z)是速度的倒数，简称为慢度，L为地震初至波的射线路径, 而T则是沿该射线路径的初至波走时.对于给定的近地表速度模型，地震初至波的射线路径L可以通过各种射线追踪方法得到，并据此射线路径计算出沿地震射线的走时.
如果给予慢度模型S(x, z)一个微小慢度扰动量ΔS(x, z)，其相应的走时变化为
|
(2) |
该方程是对走时方程(1)的一阶近似.地震初至走时层析反演则利用一阶近似方程(2)，对实际地震记录的初至时间进行反演，构建近地表的慢度模型(Zhu et al., 1992).
走时变化一阶近似方程(2)的离散形式可表达成矩阵-向量形式为
|
(3) |
式中向量ΔT为地震波初至走时的残差(实测走时减去理论走时)，矩阵A为射线路径，向量ΔS为慢度模型的修正量.向量ΔT和矩阵A为已知元素，求解线性方程组(3)可得到向量ΔS，新的慢度模型向量为S+ΔS.
一旦得到新的慢度模型S，新的射线路径既可在原有的射线路径上做相应修正，也可利用射线追踪方法重新求取.根据新的射线路径矩阵A、新的走时T和新的残差ΔT，可求得新的修正量ΔS.如此循环迭代，直到残差的最小二乘模‖ΔT‖(或者RMS数值)小于设定的门槛值为止，即获得了最终的慢度模型.
图 1中展示了层析反演静校正的基本流程.图中慢度为标量，即为速度的倒数.基本流程从上往下可分为三个区块.第一区块是建立初始模型；第二区块则是关于层析反演的循环迭代；第三区块才是依据层析反演所得到的近地表速度模型计算静校正量并进行静校正处理.
|
图 1 层析反演静校正的基本流程 Fig. 1 Flowchart of static correction by tomographic inversion |
层析反演开始之前的重要环节是建立合适的初始速度模型，该初始模型将直接影响到初至走时层析反演循环迭代的最终结果.我们根据现有的地表露头资料，结合小折射和微测井速度调查信息，构建随深度而变化的一维速度模型.虽然假设该一维速度模型适用于横向空间上的每一点，但是参照拾取的初至时间对各点的一维速度模型做了适当修正(如图 2a所示).在层析反演的循环迭代过程中，每一次模型修正之后，我们依据修正的速度模型重新进行射线追踪.重新射线追踪能够有效减小对初始模型的依赖性，但是无法改变速度模型的纵向长波长特征.因此，初始的一维速度模型需要有准确的长波长纵向速度特征，从而约束迭代反演的结果.
|
图 2 初始速度模型和反演所得的最终速度模型 (a)初始速度模型, 反映准确的速度纵向长波长特征；(b)反演所得的黄土塬近地表速度模型, 图中黑线为静校正所用的速度模型底线，黑线下方用高速替代层. Fig. 2 Initial velocity model and final velocity model from tomographic inversion (a) Initial velocity model which reflects long wavelength velocity features; (b) Loess plateau velocity model from tomographic inversion. The black line represents the bottom of the velocity model for static correction. The velocity underneath is replaced by a high velocity. |
对起伏地表的处理，我们采用的方法是所谓的等效速度充填，将速度模型向上方空中扩展，使得起伏地表成为扩展模型的一个内部地层界面.另有一种先进处理方法是采用数学变换将起伏地表拉平，而起伏地表仍然是位于模型顶部的自由表面.但是射线追踪的程函方程以及波场模拟的波动方程都需要从物理坐标变换到计算坐标中(Rao and Wang, 2013; Ma and Zhang, 2014; Rao and Wang, 2018).对于初至走时层析反演问题，采用等效速度充填扩展速度模型已经能够满足实际精度的需求.对速度模型进行网格划分时，考虑到地表起伏变化的复杂程度，应适当地将矩形单元往大面积划分，尽可能使每个网格单元中有足够多的射线穿过，从而增强反演的可靠性.图 2a中展示的初始速度模型，其横向网格宽度取15 m，纵向网格长度为8 m.
我们采用全局射线路径追踪法.该方法基于最小走时原理(Moser, 1991)，因此，射线路径追踪的实现过程是：
(1) 从源点所在节点开始，计算与其相连的所有节点的走时.
(2) 找出这些节点中走时最小的节点(该走时为沿从源点到节点的射线路径的走时，该节点射线追踪工作完成，标上“完成”标记)作为新的“源点”，走时记为ts.
(3) 计算从新的“源点”到与其相连所有节点间的走时dt，走时ts+dt较小者作为节点的新走时.
(4) 找出所有节点(不包括标有“完成”标记的点)中走时最小的节点，得到新的“完成”点与新的“源点”，走时仍然用ts表示，转到第3步.
(5) 如此不断向外追踪，直到所有接收点所在节点都标上“完成”点为止，得到一炮所有道的射线路径和走时.
3 稳定化的同步迭代重构算法线性方程组(3)式是一大型稀疏代数方程组，可以选用同步迭代重构算法(SIRT)快速求解，而不需要采用任何对矩阵广义求逆的算法.但是，它又是一个病态方程组.因此，我们给出对任意网格单元的慢度修正量的稳定化计算式为
|
(4) |
式中i指的是穿过当前矩形网格单元的第i条射线，Δti则沿该条射线的走时残差，Li是该条射线从炮点到检波点的总长度，
首先，我们需理解(4)式的物理意义.它是将每条射线的残差(Δti)加权分配到矩形慢度单元中，其权重(
其次，稳定化因子是按照整个模型中
|
(5) |
式中ε为一微小常数.该稳定化因子将有效遏制实测走时的各种误差在慢度模型重构过程中的负面影响，同时对射线稀疏区域的慢度修正量加以控制.
4 长短波长静校正方法的联合应用实际地震资料来自黄土塬覆盖区的一条二维地震测线.该地震工区的地表起伏很大，高程在920~1400 m之间.野外地震数据采集采用的是中间炸药震源激发，双边多道检波点接收的观测系统.
从实际地震记录上拾取初至波到达时间，其初至波不一定是折射波，具体是何种波型取决于局部模型特性.如果局部模型的代表性特征是均匀介质特性，最先到达地表的波是直达波.如果局部模型的代表性特征是连续渐变介质特性，初至波是回折波.只有当局部模型的代表性特征是层状介质特性时，初至波才是折射波.
初始速度模型(图 2a)展示出地表高程的剧烈变化，它将导致常规静校正方法在黄土塬覆盖区的应用遇到了困难.初至走时层析反演成功地构建黄土塬覆盖区的近地表速度模型(图 2b)，该模型清楚地反映出低降速带厚度及其速度的剧烈变化.同时我们还从模型(图 2b)上看到明显的地层分界，图中的黑线代表的是用于静校正的近地表速度模型的底线，其速度为3500 m·s-1.通常将黑线以上的速度信息用于计算静校正量.黑线的下方为高速替代层.
从层析反演的走时残差曲线(图 3a)可以看出，迭代反演在稳定快速地收敛.在第四次迭代后，已经收敛到－20 dB.迭代开始时走时残差的RMS数值为220 ms, 第四次迭代之后走时残差的RMS数值已经减小到它的10%以下，即＜22 ms.如果继续迭代的话，第9次迭代走时残差的RMS数值也仅减小到20 ms左右.由于走时拾取误差的存在，更多次迭代对走时残差的RMS数值不再会有实质意义的影响.
|
图 3 层析反演的收敛过程和射线密度分布图 (a)层析反演的收敛过程；(b)近地表速度模中的射线密度分布. Fig. 3 The convergence of the tomographic inversion and ray-path density distribution map (a) The convergence of the tomographic inversion; (b) Ray-path density distribution over the reconstructed near-surfacevelocity model. |
图 3b显示的层析模型射线密度图，色标代表的是每个网格中穿过的射线数.从图中可以看出，模型底部射线密度较低.原因是底部较深，射线进入少，但层析静校正关注的是浅层信息.可以发现，浅层中有足够多的射线穿过，由此说明反演所得近地表速度模型的可信程度.
正如图 2的速度模型所展示，黄土塬覆盖区的地表高程横向变化非常剧烈，相邻检波点的高程差异很大.由于层析反演具有统计平滑特性，根据层析反演模型计算的静校正量是长波长静校正量，应用之后仍然存在检波点的剩余静校正量.如果利用常规反射波特征进行后续的剩余静校正，不能取得满意效果.因此，需要把因巨大高程差异而导致的初至波剩余静校正量考虑进来，对计算的初至波走时进行修改.
图 4a与图 4b对比了层析反演静校正前后的炮集记录.从任意选取一对炮集记录对比可以看出，层析反演静校正之后的炮集记录(图 4b)上不仅有效反射波同相轴更加清楚，而且初至波同相轴也更加光滑和连续.而图 4c是层析静校正与初至波剩余静校正方法相结合应用之后的炮集记录.这对炮集记录上初至波剩余静校正的效果与仅仅做了层析静校正的效果相差不大，但是，在叠加剖面上将会有明显差异.
|
图 4 静校正前后的炮集记录对比 (a)层析静校正之前的炮集记录；(b)层析静校正之后炮集记录；(c)层析静校正与初至波剩余静校正方法相结合应用之后的炮集记录. Fig. 4 Comparison of shot gathers before and after tomographic static correction (a) Shot gathers before tomographic static correction; (b) Shot gathers after tomographic static correction; (c) Shot gathers obtained by tomographic static correction plus the first-arrival residual static correction. |
图 5对比了不同静校正方法处理之后的叠加剖面.首先是折射静校正处理的结果(图 5a)，与层析静校正处理的结果(图 5b)的对比.层析静校正处理的深层同相轴连续性好，反射波组特征清晰.这是因为层析反演构建了较为准确的近地表速度模型.
|
图 5 黄土塬覆盖区不同静校正方法的叠加剖面对比 (a)常规折射静校正后叠加剖面; (b)层析静校正后叠加剖面; (c)层析静校正与初至波剩余静校正方法相结合应用之后的叠加剖面. Fig. 5 Comparison of stack profiles by various static corrections in the loess plateau (a) Conventional refraction static correction; (b) Tomographic static correction; (c) Tomographic static correction combined with first-arrival residual statics correction. |
图 5c是综合应用了初至波层析静校正和剩余静校正之后得到的叠加剖面.由于剩余静校正考虑了高程的剧烈变化，叠加剖面质量明显提升.不但目的层反射的信噪比增高，构造形态更加清楚，而且最终剖面的浅层地震信息更加丰富.
5 结论黄土塬表层覆盖巨厚黄土，地表起伏高差较大，地震静校正问题严重.而且黄土塬覆盖区的潜水面较深，常规静校正方法无法取得理想效果.本文针对黄土塬覆盖区的静校正难题，研究层析反演静校正方法在黄土塬地区的适用性和可靠性.层析反演静校正利用地震波初至走时数据、通过迭代反演的方法构建速度模型，本文的迭代反演修正量的计算采用同步迭代重构算法.走时方程组是病态线性方程组，因此本文对同步迭代重构算法进行了改进，使得层析反演的迭代过程趋于稳定.
层析静校正具有统计平滑特性，可是黄土塬地区地表高程的横向变化剧烈，相邻检波点的高差及其静校正量有时差异很大.因此，在运用层析静校正求取长波长静校正量的同时，还需计算初至波的剩余静校正量.实例证明，综合应用层析静校正和剩余静校正方法，同时考虑了黄土塬覆盖区的长波长和短波长的静校正量，能够有效地解决实际地震资料的静校正问题.
Abbas A, Sarosh B, Shah S F. 2018. Computation of residual statics in complex geological areas using simulated annealing method based on minimization of objective function. Journal of Geophysics and Engineering, 15(6): 2577-2585. DOI:10.1088/1742-2140/aadb21 |
Cha W F, Yu X L. 2014. The optimal selection and comprehensive application of the static correction in the complex shallow surface structure areas. Progress in Geophysics (in Chinese), 29(2): 958-965. DOI:10.6038/pg20140265 |
Cox M. 1999. Static Corrections for Seismic Reflection Surveys. Tulsa: Society of Exploration Geophysicists.
|
Dahl-Jensen T. 1989. Static corrections on crystalline rock. Geophysical Prospecting, 37(5): 467-478. DOI:10.1111/j.1365-2478.1989.tb02218.x |
Jing X L, Yang C C, Li Y M, et al. 2002. A global optimized algorithm for seismic residual statics corrections. Chinese Journal of Geophysics (in Chinese), 45(5): 707-713. |
Le Meur D, Poulain G, Bertin F, et al. 2011. Monte-Carlo statics on P-P or P-Sv Wide-Azimuth data.//81st Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
|
Lin Y H. 2003. Hybrid optimization of static estimation in complex topography. Chinese Journal of Geophysics (in Chinese), 46(1): 101-106. |
Ma T, Zhang Z. 2014. Calculating ray paths for first-arrival travel times using a topography-dependent eikonal equation solver. Bulletin of the Seismological Society of America, 104(3): 1501-1517. DOI:10.1785/0120130172 |
Marsden D. 1993. Static correction:a review, Part 1. The Leading Edge, 12(1): 43-49. DOI:10.1190/1.1436912 |
Moser T J. 1991. Shortest path calculation of seismic rays. Geophysics, 56(1): 59-67. DOI:10.1190/1.1442958 |
Rao Y, Wang Y H. 2013. Seismic waveform simulation with pseudo-orthogonal grids for irregular topographic models. Geophysical Journal International, 194(3): 1778-1788. DOI:10.1093/gji/ggt190 |
Rao Y, Wang Y H. 2018. Seismic waveform simulation for models with fluctuating interfaces. Scientific Reports, 8: 3098. DOI:10.1038/s41598-018-20992-z |
Rothman D H. 1985. Nonlinear inversion, statistical mechanics, and residual statics estimation. Geophysics, 50(12): 2784-2796. DOI:10.1190/1.1441899 |
Rothman D H. 1986. Automatic estimation of large residual statics corrections. Geophysics, 51(2): 332-346. DOI:10.1190/1.1442092 |
Wiggins R A, Larner K L, Wisecup R D. 1976. Residual statics analysis as a general linear inverse problem. Geophysics, 41(5): 922-938. DOI:10.1190/1.1440672 |
Yang C Z, Jin D M, Liang D W, et al. 2016. Method of identification and solution about long-wavelength static correction:a case study from the seismic data in Ordos loess plateau. Progress in Geophysics (in Chinese), 31(5): 2212-2218. DOI:10.6038/pg20160545 |
Yang K, Cheng J B, Liu Y Z, et al. 2007. A study on the application of the 3-D wave equation datuming. Chinese Journal of Geophysics (in Chinese), 50(4): 1232-1240. |
Zhu X H, Sixta D P, Angstman B G. 1992. Tomostatics:turning-ray tomography+static corrections. The Leading Edge, 11(12): 15-23. DOI:10.1190/1.1436864 |
Zhang L, Yang Q Y, Zhang B, et al. 2017. Tomography inversion by first breaks in areas with complex near surface. Progress in Geophysics (in Chinese), 32(2): 816-821. DOI:10.6038/pg20170249 |
查文锋, 于小磊. 2014. 浅表层结构复杂区静校正的优化选取及应用. 地球物理学进展, 29(2): 958-965. DOI:10.6038/pg20140265 |
井西利, 杨长春, 李幼铭, 等. 2002. 地震静校正全局最优化问题的求解. 地球物理学报, 45(5): 707-713. DOI:10.3321/j.issn:0001-5733.2002.05.013 |
林依华. 2003. 复杂地形条件下静校正的综合寻优. 地球物理学报, 46(1): 101-106. DOI:10.3321/j.issn:0001-5733.2003.01.016 |
杨城增, 金东民, 梁殿文, 等. 2016. 长波长静校正问题的识别与解决方法——以鄂尔多斯黄土塬地震资料为例. 地球物理学进展, 31(5): 2212-2218. DOI:10.6038/pg20160545 |
杨锴, 程玖兵, 刘玉柱, 等. 2007. 三维波动方程基准面校正方法的应用研究. 地球物理学报, 50(4): 1232-1240. DOI:10.3321/j.issn:0001-5733.2007.04.033 |
张林, 杨勤勇, 张兵, 等. 2017. 复杂近地表初至波层析反演静校正技术研究. 地球物理学进展, 32(2): 816-821. DOI:10.6038/pg20170249 |