CN109633749B - 基于散射积分法的非线性菲涅尔体地震走时层析成像方法 - Google Patents

Abstract

本发明涉及一种基于散射积分法的非线性菲涅尔体地震走时层析成像方法，包括以下步骤：1)对原始地震数据进行去噪和反褶积的预处理；2)在预处理过的地震数据炮记录上进行初至拾取，获得拾取走时；3)根据地下先验信息建立初始表层速度模型，并设定反演参数，包括反演频率范围和间隔；4)根据初至波到达时和初始表层速度模型进行散射积分法非线性菲涅尔体地震层析成像反演，获得最终表层速度模型并成像。与现有技术相比，本发明具有占用的内存小、计算效率高、易于并行、方便预条件、反演精度高、过程稳定等优点。

Description

基于散射积分法的非线性菲涅尔体地震走时层析成像方法

技术领域

本发明涉及地震走时层析成像领域，尤其是涉及一种基于散射积分法的非线 性菲涅尔体地震走时层析成像方法。

背景技术

射线走时层析被广泛应用在天然地震学和勘探地震学中(Sheriff&Geldart,1982；Dziewonski,1984；Nolet,1987；Pulliam et al.,1993；Billette&Lambaré,1998)。但是由于该方法基于高频近似，所以其反演分辨率的提高受到了限制。因此， Yomogida(1992)提出了与频率有关的层析方法，Vasco and Majer(1993)将波路径走 时敏感核函数应用到井间层析反演中并获得了比传统射线走时层析具有更高分辨 率的反演结果，Marquering等(1998,1999)和Dahlen等(2000)进一步提出了有限频 层析并深入研究了其对应核函数。随后，有限频层析被广泛应用到区域的和全球的 地震数据中(Montelli etal.,2004；Yoshizawa&Kennett,2004；Zhou et al.,2006； Sigloch et al.,2008；Tian etal.,2009)。为了降低内存的占用量，Spetzler和Snieder (2004)将敏感核函数限制在第一菲涅尔体中并提出了菲涅尔体层析，Liu(2009)在此 基础上进一步分析了带限敏感核函数的特征。菲涅尔体层析方法在近地表速度反演 和二氧化碳储藏的时移地震监测中的应用展现出了其巨大的潜力。但是传统的FVT 需要求解大规模病态矩阵，所以很占内存并且常常不稳定。Tromp等(2005)和 Taillandier等(2009)利用伴随状态法构建梯度并分别将其应用于有限频层析和射线 走时层析中，然而伴随状态法很难实施预条件。散射积分法(SI)(Chen et al.,2007) 则是另一种计算梯度的办法。该方法通过显式地计算核函数，并与走时残差向量相 乘实现梯度的计算。它的计算效率取决于观测系统的设置，当炮点数多于检波点数 时，它相比于伴随状态法具有计算量上的优势(Chen,2007；Liu et al.,2015)。但与传 统的层析方法一样，散射积分法面临内存占用大的问题，特别是在利用Hessian矩 阵时问题更加突出。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于散射积 分法的非线性菲涅尔体地震走时层析成像方法。

本发明的目的可以通过以下技术方案来实现：

一种基于散射积分法的非线性菲涅尔体地震走时层析成像方法，包括以下步骤：

1)对原始地震数据进行去噪和反褶积的预处理；

2)在预处理过的地震数据炮记录上进行初至拾取，获得拾取走时；

3)根据地下先验信息建立初始表层速度模型，并设定反演参数，包括反演频 率范围和间隔；

4)根据初至波到达时和初始表层速度模型进行散射积分法非线性菲涅尔体地 震层析成像反演，获得最终表层速度模型并成像。

所述的步骤4)具体包括以下步骤：

41)在当前表层速度模型上进行射线追踪获取理论合成初至走时以及模型空间的旅行时场；

42)判断理论合成走时与拾取走时是否匹配，若是，则将当前表层速度模型作 为最终表层速度模型，若否，则执行步骤43)；

通过二范数形式的目标函数是否小于实验最初设定的期望值，如果大于期望值则认为不满足匹配要求，执行步骤43)，如果小于预设的期望值，则认为满足要 求，停止迭代；

43)采用LU分解算法获取所有频率的炮点端波场与检波点端格林函数；

44)根据炮点端波场和格林函数计算每一个炮检对所对应的核函数，并根据旅 行时场确定菲涅尔体的范围，获得带核函数特征的菲涅尔体；

45)通过核函数-标量乘获得每一个炮检对所对应的梯度，累加形成整个观测 系统所对应的梯度，同时计算核函数自身元素的模平方并沿着所有炮检对累加求和 获得梯度的预条件算子；

46)计算预条件最速下降方向、方向更新步长和模型更新量后更新当前表层速 度模型；

47)重复步骤41-46)，直到获得匹配的最终表层速度模型。

所述的步骤41)具体包括以下步骤：

411)在当前表层速度模型网格剖分的网格边界上设定多个关键点；

412)自激发点开始以矩形为波前，按惠更斯原理逐层向外扩展，直至扫描完 整个模型空间，再由外向里收缩，直至扫描到激发点，如此重复直至模型空间的最 小旅行时场不变为止，在扩展或收缩过程中，次级源来自于两个关键点之间，两个 关键点之间的走时通过线性插值得到；

413)检波点处的初至波到达时即为理论合成初至走时。

所述的步骤44)具体包括以下步骤：

441)计算菲涅尔体核函数K(r,ω|g,s)，计算式为：

其中，ω为圆频率、Im表示取虚部、G₀(g,r)为在背景介质v₀(r)中空间位置r 处激发g点接收到的格林函数，u₀(r,s)为空间位置s处激发r点接收到的入射波场， u₀(g,s)为空间位置s处激发g点接收到的入射波场；

412)根据公式确定菲涅尔体的范围，其中τ(r,s)表示激发点s至空间点r的最小走时，τ(r,g)表示接收点g至空间点r的最小走时， τ(g,s)表示激发点s至接收点g的最小走时，T表示主频的周期。

所述的步骤45)中，梯度与预条件算子的计算步骤如下：

451)获取每一个炮检对所对应的菲涅尔体以及菲涅尔体与此炮检对所对应的 走时残差乘，得到此炮检对所对应的梯度以及菲涅尔体向量的元素平方向量；

452)将所有的炮检对梯度累加求和得到全局梯度K^TΔt，具体计算式为：

其中，k_ij为每一对炮检对对应的Freach核函数，Δt_i为旅行时差；

453)将所有的炮检对的核函数元素取模平方，并进行累加，得到全局梯度的 预条件算子H₀，具体计算式为：

所述的步骤46)中，预条件最速下降方向p的计算式为：

所述的步骤46)中，方向更新步长的计算方法为：

根据先验信息获取初始模型与真实地下模型的速度差，将速度差除以预先设定的迭代次数，得到每次迭代可更新的最大速度值Δv_max，利用待定系数法求取步长 t，使得max{|p·t|}＝Δv_max。

与现有技术相比，本发明具有以下优点：

一、占用的内存小：该方法无需存储核函数矩阵即可方便地实现目标函数梯度 的计算，在任何时刻只需存储单个菲涅尔体核函数；

二、计算效率高，易于并行：该方法将大规模的核函数-向量乘运算表示为具 有明确物理含义的向量-标量乘的累加运算，无需利用SVD或者LSQR方法解大规模 矩阵方程组，计算时间大大减少；

三、方便预条件，使得计算效率进一步提高，反演精度也大大提高：该方法无 需存储Hessian矩阵即可方便实现预条件，预条件算子的应用可以加快收敛速度， 并能够获得地下深部介质的信息；

四、反演过程更稳定：优化的梯度导引算法和预条件算子的使用使得反演过程 更稳定。

附图说明

图1为本发明的流程图。

图2为本发明的硬件结构示意图。

图3为实施例1的真实理论模型图。

图4为实施例1的常梯度初始模型图。

图5为实施例1的真实理论模型上半图。

图6为实施例1的地震记录垂直分量图，其中，图(6a)为地表水平9km的 地震记录垂直分量图，图(6b)为地表水平26km的地震记录垂直分量图。

图7为实施例1的本发明(SI-FVT)反演结果图。

图8为实施例1的SIRT-FVT反演结果图。

图9为实施例1的真实模型(黑色实线)，初始模型(线点)与SI-FVT(线点 点)和SIRT-FVT(间隔线段)在地下不同深度处的速度对比图，其中，图(9a) 为在地下100m深度处的速度对比图，图(9b)为在地下200m深度处的速度对比 图。

图10为实施例1的真实模型，SI-FVT和SIRT-FVT得到的反演结果对应初至 走时的对比图，可以看到三者几乎重叠在一起。

图11为实施例1的SI-FVT(虚点)和SIRT-FVT(线线点)收敛曲线图。

图12为实施例2的初始速度模型图。

图13为实施例2的本发明(SI-FVT)反演结果图。

图14为实施例2的SIRT-FVT反演结果图。

图15为实施例2的SI-FVT和SIRT-FVT在不同水平位置的垂直速度对比图， 其中，图(15a)为在3km水平位置的垂直速度对比图，图(15b)为在6km水平 位置的垂直速度对比图，图(15c)为在9km水平位置的垂直速度对比图。

图16为实施例2的真实拾取，SI-FVT和SIRT-FVT得到的反演结果对应初至 走时的对比图，也可以看到几乎重叠在一起。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

在我国西部，地震勘探逐渐从前期的沙漠、戈壁滩地区扩展到地表条件更为复 杂的山地、山前带，如塔西南、库车、吐哈、青海、陕甘宁黄土塬等。山地山前地 区地形起伏剧烈、表层结构复杂、速度横向变化大，折射界面不稳定或者不存在， 反映在地震资料上表现为严重的噪声干扰问题、能量失衡问题、静校正问题等。其 中，静校正问题是首要问题，这一问题的解决是解决其它问题的前提。常规的静校 正技术很难适用于复杂地表地区，因此，准确估算表层速度模型，借此消除复杂表 层因素在地震资料上产生的静校正问题，是野外勘探和室内资料处理的共同任务。

在复杂地表地区，地形起伏和近地表速度强烈变化会给常规采集和处理技术带来许多问题，然而，通过使用初至波层析成像技术能够较好地回避或解决这些问题。 通过使用本发明提出的基于改进的散射积分法的非线性菲涅尔体走时层析，相比于 常规初至波走时层析更具有优势。

实施例1：

本实施例将以二维复杂起伏地表理论模型作为真实模型(如图3所示)，该模 型共有4001×151个网格，网格间距为10m×10m，速度最大最小值分别为800m/s 和4300m/s。在该模型上进行弹性波正演模拟，共正演751炮，炮间距为40m，第 一炮在5000m的位置，每炮有301个检波器，分布在炮点两侧，道间距为20m。 因此该正演记录的最大偏移距为3000m，最小偏移距为0m。在应用本发明提出的 基于改进的散射积分法的非线性菲涅尔体层析(SI-FVT)的同时，基于SIRT法的 非线性菲涅尔体层析(SIRT-FVT)也同时被应用，以对比突出本发明的有效性和 优越性。由于初至波走时层析只能得到近地表速度结构，所以对比结果只展示速度 模型的上半部分，如图5所示。

具体实施方式如下：

(1)数据采集器1采集地震波信号，并对原始地震数据进行去噪、反褶积等 非破坏走时的预处理(如图6所示)；

(2)数据采集器1将步骤(1)预处理后的地震波数据逐道输入到处理器中， 在处理过的地震数据炮记录上进行初至拾取，以获得初至到达时；

(3)输入设备4根据地下先验信息建立初始速度模型(如图4所示)，并设 定反演频率范围与间隔等反演参数；

(4)处理器2在当前模型上进行射线追踪以获得理论合成初至走时及模型空 间的旅行时场；

(5)处理器2判断理论合成走时与拾取走时的匹配程度，如果“吻合”则保 留当前模型并退出，并通过显示器显示模型(如图7、8、9、10所示)，否则继续 执行以下步骤；

(6)处理器2计算并存储所有频率的炮点端波场与检波点端格林函数；

(7)处理器2根据步骤(6)的炮点端波场和格林函数计算每一个炮检对所对 应的核函数，并利用步骤(4)输出的旅行时场圈定菲涅尔体范围，以获得带核函 数特征的菲涅尔体(以下简称菲涅尔体)；

(8)处理器2通过核函数-标量乘获得每一个炮检对所对应的梯度，并累加形 成整个观测系统所对应的梯度，同时计算核函数自身元素的模平方并沿着所有炮检 对累加求和获得梯度的预条件算子；

(9)处理器2计算预条件最速下降方向；

(10)处理器2计算方向更新步长；

(11)处理器2计算模型更新量，并更新模型；

(12)重复步骤(4)-(11)；

利用SI-FVT和SIRT-FVT进行反演得到的结果分别如图7和图8所示。由于 缺少观测记录覆盖，所以模型两侧的反演结果不如中间部分准确。对真实模型，初 始模型以及两种方法所得到的反演结果在不同深度进行切片对比，如图9所示。不 同的反演结果对应的初至波走时如图10所示。可以看到两种反演方法得到的反演 结果所预测的初至走时与真实模型的初至走时有很好的吻合度。但从图9可以看出， SI-FVT的反演结果相比于SIRT-FVT有着较高的分辨率和精度。这可能是因为 SI-FVT稳定性较好，并且利用了预条件算子。在物理上，预条件算子能够补偿照 明，从而反演更深部的信息并且提高反演的精确度。图11为不同算法的收敛曲线。 从图上可以看出SI-FVT有着更好的稳定性，并且在最后一步之后SI-FVT的目标 函数值更小。表1展示了两种算法每轮循环的内存占用量以及计算时间，显然 SI-FVT的计算效率要高于SIRT-FVT，而其内存占用量要比SIRT-FVT少一个数量 级，体现了本发明的优势。

表1：两种方法每轮循环所需内存量和计算时间对比

实施例2：

本实施例将本发明提出的基于改进的散射积分法的非线性菲涅尔体走时层析 应用到我国西部四川盆地获取的实际资料中。模型共有2957×144个网格，网格间 距15m×15m。共500炮，炮间距60m，第一炮在地表水平位置7185m处。每炮 480道，均匀分布在炮点两端，道间距30m。该观测系统最大偏移距7185m，最 小偏移距15m。具体实施流程与实施例1相似，其中步骤(3)中输入的初始模型 如图12所示，显示器输出的结果如图13～图16所示。

分析SI-FVT和SIRT-RT(基于SIRT法的射线走时层析)的反演结果(图13、 图14)，SI-FVT的反演结果浅部分辨率更高，尤其在箭头所指的位置。其次，最 左边圈出的位置二者的反演结果明显不同，SI-FVT的反演结果更能满足地表一致 性的先验信息，因此更加合理。从方框标出的位置也可以明显看到SIRT-RT的反 演结果存在“脚印”，图15的箭头也标示出了这些“脚印”的存在。图16展示了 不同的反演结果对应的初至波走时，可以看到SI-FVT对应的初至波走时更贴近真 实拾取的走时。以上证明了FVT相比于RT能获得更高分辨率的结果，尤其在深 部位置。

Claims (6)

1.一种基于散射积分法的非线性菲涅尔体地震走时层析成像方法，其特征在于，包括以下步骤：

1)对原始地震数据进行去噪和反褶积的预处理；

2)在预处理过的地震数据炮记录上进行初至拾取，获得拾取走时；

3)根据地下先验信息建立初始表层速度模型，并设定反演参数，包括反演频率范围和间隔；

4)根据初至波到达时和初始表层速度模型进行散射积分法非线性菲涅尔体地震层析成像反演，获得最终表层速度模型并成像，具体包括以下步骤：

41)在当前表层速度模型上进行射线追踪获取理论合成初至走时以及模型空间的旅行时场；

42)判断理论合成走时与拾取走时是否匹配，若是，则将当前表层速度模型作为最终表层速度模型，若否，则执行步骤43)；

43)采用LU分解算法获取所有频率的炮点端波场与检波点端格林函数；

44)根据炮点端波场和格林函数计算每一个炮检对所对应的核函数，并根据旅行时场确定菲涅尔体的范围，获得带核函数特征的菲涅尔体；

45)通过核函数-标量乘获得每一个炮检对所对应的梯度，累加形成整个观测系统所对应的梯度，同时计算核函数自身元素的模平方并沿着所有炮检对累加求和获得梯度的预条件算子；

46)计算预条件最速下降方向、方向更新步长和模型更新量后更新当前表层速度模型；

47)重复步骤41-46)，直到获得匹配的最终表层速度模型。

2.根据权利要求1所述的一种基于散射积分法的非线性菲涅尔体地震走时层析成像方法，其特征在于，所述的步骤41)具体包括以下步骤：

411)在当前表层速度模型网格剖分的网格边界上设定多个关键点；

412)自激发点开始以矩形为波前，按惠更斯原理逐层向外扩展，直至扫描完整个模型空间，再由外向里收缩，直至扫描到激发点，如此重复直至模型空间的最小旅行时场不变为止，在扩展或收缩过程中，次级源来自于两个关键点之间，两个关键点之间的走时通过线性插值得到；

413)检波点处的初至波到达时即为理论合成初至走时。

3.根据权利要求2所述的一种基于散射积分法的非线性菲涅尔体地震走时层析成像方法，其特征在于，所述的步骤44)具体包括以下步骤：

441)计算菲涅尔体核函数K(r，ω|g，s)，计算式为：

其中，ω为圆频率、Im表示取虚部、G₀(g，r)为在背景介质v₀(r)中空间位置r处激发g点接收到的格林函数，u₀(r，s)为空间位置s处激发r点接收到的入射波场，u₀(g，s)为空间位置s处激发g点接收到的入射波场；

412)根据公式

确定菲涅尔体的范围，其中τ(r，s)表示激发点s至空间点r的最小走时，τ(r，g)表示接收点g至空间点r的最小走时，τ(g，s)表示激发点s至接收点g的最小走时，T表示主频的周期。

4.根据权利要求2所述的一种基于散射积分法的非线性菲涅尔体地震走时层析成像方法，其特征在于，所述的步骤45)中，梯度与预条件算子的计算步骤如下：

451)获取每一个炮检对所对应的菲涅尔体以及菲涅尔体与此炮检对所对应的走时残差乘，得到此炮检对所对应的梯度以及菲涅尔体向量的元素平方向量；

452)将所有的炮检对梯度累加求和得到全局梯度K^TΔt，具体计算式为：

其中，k_ij为每一对炮检对对应的Freach核函数，Δt_i为旅行时差；

453)将所有的炮检对的核函数元素取模平方，并进行累加，得到全局梯度的预条件算子H₀，具体计算式为：

5.根据权利要求4所述的一种基于散射积分法的非线性菲涅尔体地震走时层析成像方法，其特征在于，所述的步骤46)中，预条件最速下降方向p的计算式为：

6.根据权利要求5所述的一种基于散射积分法的非线性菲涅尔体地震走时层析成像方法，其特征在于，所述的步骤46)中，方向更新步长的计算方法为：

根据先验信息获取初始模型与真实地下模型的速度差，将速度差除以预先设定的迭代次数，得到每次迭代可更新的最大速度值Δv_max，利用待定系数法求取步长t，使得max{|p·t|}＝Δv_max。

Images