全波形反演中最优化传输函数研究

胡光辉,贺伟光

(中石化石油物探技术研究院有限公司,江苏南京211103)

全波形反演(FWI)将叠前地震数据信息映射为地下结构或介质参数分布,是目前地震学中最强大的建模成像工具之一。在勘探地球物理学中,早在1975年就建立了从地面地震数据推断内部结构的概念[1-2],但FWI的理论建立者仍被认为是LAILLY[3]和TARANTOLA[4]。对于全波形反演一词的理解存在两个观点:第1个观点坚持认为只有利用整个地震炮集数据的波形信息才能称为全波形反演(full指的是整个地震炮集数据);第2个观点更宽泛一些,把所选择的地震事件信息完全利用即可成为全波形反演(full暗含的意思是fully)。无论是哪种观点,都认为波的传播需要经过求解波动方程来模拟,利用其它退化的模拟方法(比如射线方程法、单程波方程法)都不能称为全波形反演。在全波形反演的早期发展中,由于计算能力有限以及缺乏低频和长偏移距观测数据,因而其发展十分缓慢。自2000年以来,由于全波形反演理论上的发展和计算能力的进步,该方法逐步发展成为现代工业生产中的常规地震资料处理方法之一,尤其是海上地震资料。然而,无论是在全波形反演的早期发展阶段还是如今的工业应用阶段,周跳(cycle skipping)问题依然没有完全解决[5-6]。尽管有一些精巧的反演策略[7](频率逐渐增加和偏移距逐渐增加),但是传统的最小二乘函数(L2函数)在许多情况下依然会陷入局部极值。

从数学上讲,全波形反演是在弹性动力学方程约束下的局部优化问题。数学上的处理有两种:第1种思路是构建拉格朗日函数;第2种思路是建立增广拉格朗日函数。增广拉格朗日函数是将弹性动力学方程作为二次项来构建优化函数[8]。从数学界来看,增广拉格朗日函数是凸性良好的优化函数。波动方程以二次型的方式引入增广拉格朗日函数并不要求波动方程被精确地满足。这种处理方式相当于增大了波形反演的解的搜索空间,同时增强了方程的凸性。但是,增广拉格朗日函数的计算要求很高,原因是涉及模拟算子和伴随算子的乘积。在频率域,计算代价可以接受,但在时间域则相当困难。因此,增广拉格朗日函数在工业生产中很少使用。全波形反演的主流框架一直是第1个思路,PLESSIX[9]对该框架进行了综述。在拉格朗日框架中,作为约束的弹性动力学方程被线性地引入到优化函数中,之后可以方便地得到伴随场震源和梯度的计算。拉格朗日框架要求在每次迭代时都能准确地求解弹性动力学方程。等价地说,在全波形反演的整个迭代过程中,约束方程(弹性动力学方程)始终得到满足。这限制了未知参数的搜索空间,因此导致拉格朗日框架的凸性不如增广拉格朗日框架。但是,拉格朗日框架因为计算代价较低而成为学术界和工业界的主流方式。在拉格朗日函数框架方案中,避免周跳的主要方法是设计凸目标函数,如L1范数[10]、混合或Huber函数[11-13]、零延迟互相关[14-17](相当于归一化的L2函数)、包络函数[18-24]、惩罚互相关函数[25]、基于反褶积的函数[26]和最优传输函数[27-33]。其中,最优化传输函数的优势是其比较的是全局地震数据,而不是逐点比较地震数据。这种特性保持了地震数据的事件间和检波点端记录数据的相干性,从而使得较传统L2函数能更稳定地反演。本研究采用拉格朗日函数的框架,综述了几种最优传输函数在全波形反演中的应用。

最优化传输起源于法国的一个工程问题,即如何将沙子从源头位置运输出来,以最少的能耗来填补施工区域的洞。这个问题被MONGE[34]用数学方法表述为一个L1最小化问题。在一些研究中,它可能被命名为EMD(Earth Mover Distance)。方程的答案可以通过求解非线性方程,这种非线性方程也称为Monge-Ampere方程。将1范数推广到任何其它范数称为p-Wasserstein函数,其中,p表示范数。证明Monge方程解的存在性,以及求解都是很困难的。KANTOROVICH[35]取得了突破,将Monge公式扩展到更高维空间会得到线性问题;具体来说,这是线性分配问题。FERRADANS等[36]和LELLMANN等[37]在图像处理中引入最优化传输函数,该函数是一种自动且正确的匹配模式。最优化传输函数目前已广泛应用于机器学习或人工智能(AI)等领域。最优传输问题根据分布是连续的还是离散的[38]分为3类,即连续(两个分布都是连续的)、半离散(一个是连续的,另一个是离散的)和离散的(两者都是离散的)。连续问题可以用Monge-Ampere方程或流体动力学方程来解决。半离散问题利用计算几何和牛顿优化技术可以快速解决。离散问题用正则化最优传输求解器或标准线性赋值求解器求解。

直接计算最优传输问题的运算量巨大。例如,当将n个点移动到n个目标位置时,每个点都可以移动到任何位置。因此,未知变量的数目为n2。KANTOROVICH[35]将原始的最优化传输问题转换到对偶空间中,未知量的个数变为n,极大地减少了离散情形下未知变量的个数。在对偶空间中,这可以用标准的线性方程求解器来方便地求解。同时,也有一些专门针对线性分配的算法加速收敛,比如匈牙利算法[39]或拍卖算法[40]。在图像比较领域,在传输代价函数中添加关于传输算子的正则化被广泛采用[41]。为了提高计算效率,人们发展了更为精巧和复杂的计算技术,如多层次采样策略、多尺度正则化策略。在多层次采样策略中,将整个数据集采样到多个层次。将原始数据集定义为零级,沿每个维度以2的比率进行降采样得到第1级(对于2D问题,这个数字现在是原始数据集的0.25倍)。第1级再次重新采样以获得下一级,最后一级是最小的数据集合。数据分级采样之后,求解器从最小的数据集开始(也是最终的数据集)。最后一级的传输平面计算出来之后,进行插值,成为前一级的初始会传输平面。这个过程一直持续到零级。在每个层面上,都采用多尺度的正则化策略来计算传输平面。在各级计算中,正则化方法相应地从相对较大的值开始,逐渐减小到接近于零。在传输平面代表的映射中,每个点会被分割到几个位置上。在零级数据集上,传输平面上会有一些非常集中的能量分布,每个能量团的中心对应于一个一一映射的解。正则化的引入使得无法得到单点对单点的映射,但是可以定义为质心坐标强制传输平面成为一对一的映射。

最优化传输并不能直接应用到地震数据处理中。地震信号是上、下震荡,有正、有负的。而最优化传输处理的是正信号,且两个比较对象的质量要相等。因此需要对数据进行预处理,使得地震信号满足非负且质量相等的条件。依据数据预处理的方式,目前有3种最优化传输公式在全波形反演中得到了很好的应用,即2-Wasserstein函数[27-29],KR-OT函数[30-32]和图空间的OT函数GSOT[33]。2-Wasserstein函数是基于范数2准则的最优传输函数的公式。它是第一个被引入到地球物理学的最优传输函数[27](W22函数,第一个2意味着2范数,第二个2意味着平方)。它是一维的算法,每次对比一对地震道信号。每一个合成地震道和观测地震道都应该经过预处理为正的并归一化(保证质量相等)。最简单的预处理是减去最小值并除以其积分。该公式的独特优点是存在解析解,计算速度快且稳定。扩展到其它规范(p-Wasserstein函数)是很直接的。缺点是预处理对理论凸性略有破坏。在极端情形下,可以假设地震信号减去负无穷大,在此情形下,最优化传输函数退化为常规的最小二乘误差。

第2个引入地球物理学的最优化传输是由METIVIER等[30]提出的KR-OT。它是基于1范数(1-Wasserstein函数)的最优化传输函数,并且是在对偶空间中表示。KR-OT在一维、二维和三维空间中都具有相似的数学结构。二维全波形反演的案例表明,2D KR-OT明显优于1D KR-OT。正如METIVIER所解释的那样,1D KR-OT反演的退化是忽略了接收器间的相干性。在KR-OT的公式中,最优化传输函数只依赖于两个数据残差(地震数据可以直接发送到求解器,而不需要移动或归一化),这是优点也是缺点。弊端在于所有的地震事件从一开始就混合在一起,当初始合成数据与观测数据相差较远时,映射可能无法正确匹配地震事件。KR-OT的局限之一是采样在时间轴和空间轴上应该是均匀的。至于将KR-OT扩展到其它领域,如超声成像[42],原则上是不合适的,因为接收器并不是均匀分布的。

上述两个最优化传输公式都属于连续传输问题。勘探地球物理学中的第3种最优传输函数是图空间最优传输函数GS-OT,这是离散传输问题。将一维离散信号扩展为时间-振幅二维图中的点集[33,43-44]。这种处理方法自然满足了数据非负的要求。同时,如果用相同数量的点离散地震数据,质量守恒也是自然满足的。从原理上将沿着任何一个维度都可以任意采样,去除了均匀采样的要求。GS-OT仅取决于点的数量,即在一个地震道中的离散采样点。所以直接扩展到高维空间必然面对计算量巨大的挑战。

最近,XU等[45]开发了基于图空间公式的高维最优传输函数。虽然这项研究是在全波形反演中的应用,但作者最初的研究动机是要匹配两个彩色点云。每个颜色点都与一个向量(R,G,B)相关联,但最优化传输只针对标量函数。从数学上讲,XU等[45]开发的是比较两个向量场的最优化传输函数。该研究进行了二维和三维实验。正如作者解释的那样,当点的数量少于1000时,计算非常快;但对于超过10000个点的彩色点云来说,计算时间过长。因为该研究将整个地震图分成互相不重叠的小区域(二维地震图是长方形,三维是长方体),大幅度降低了计算成本。依据在于,地震图的对比中,合成数据的地震事件是不应当与观测事件相差太远。所以划成小区域后,数据匹配依然存在合理性。但是在图像对比中,整张图可能会旋转或者倒转,在这种情形下不能划分小区域。经过数据区域划分处理后,将高维最优传输函数纳入全波形反演中变得可行和现实。该研究还将熵变换应用于空间坐标和时间坐标。虽然在一些模型中是否采用熵变换并不影响反演,但有趣的是,在Chevron 2014盲测中,应用熵变换的反演结果优于没有熵变换的反演。

最优化传输函数是在原始空间或对偶空间中求解。在数学领域,一种新的基于混合算法的最优传输方法被提出。该方法是基于2-Wasserstein函数的二维算法,不再是一维单道处理。这是由JACOBS等[46]开发的一种新颖的方法来解决最优化传输函数。首先,将原始形式转换到对偶空间,在数学上就得到了两个等价的方程:一个在原始空间,另一个在对偶空间。新颖的部分在于第2步。JACOBS等[46]只迭代一步原始解;然后跳到对偶空间迭代一步;之后再回到原始空间继续循环,直到得到一个收敛的结果。JACOBS等[46]将其命名为往复法(back and forth method,BFM)。

最初的最小二乘目标函数是逐点比较的算法,后来兴起的从互相关到最优传输的目标函数更多是基于全局比较算法。因此,可以更好地利用相邻点之间甚至地震事件之间的相干性。这很好地提高了目标函数的凸性和地下结构的一致性。本研究旨在提供对全波形反演中最优化传输函数的介绍,解释了目前引入全波形反演的几种最优化传输公式,包括W2(1D版本,2-Wasserstein函数)、BFM、KR-OT、GS-1D和GS-2D(高维最优传输函数)。

1 方法理论

本文重点在于讨论最优化传输函数,但为了说明与全波形反演的结合应用,先简单介绍全波形反演的实现。之后重点介绍5个最优化传输函数的基本原理。W2是1D 2-Wasserstein函数,BFM是2D或者3D的2-Wasserstein函数,KR-OT是2D或者3D的1-Wasserstein函数,GS-1D是一维信号的图空间最优化传输,GS-2D是多维信号的图空间最优化传输。

1.1 全波形反演

FWI已经发展了几十年,现在是地震资料处理业界的常规技术。FWI包含一个目标函数和一个约束方程(又称为控制方程)。目标函数是刻画合成地震数据与观测地震数据直接偏差的函数,控制方程是描述声波或弹性地震波传播的方程。最早也是最常见的目标函数是最小二乘函数(L2函数),即:

(1)

式中:p(t)是合成数据;d(t)是观测数据,求和针对所有的震源和接收器。目标函数中应用了加权函数W(t),因为并不是所有的地震事件都被利用。例如,在应用声波方程求解器处理陆地地震数据时,需要预先去除观测数据中的面波和横波。目标函数的优化通常是基于局部优化求解器,模型更新的搜索方向是基于梯度的。梯度的计算通常是采用伴随方[9]。正传波场与伴随波场之间的累积求和为梯度。计算梯度只需要花费2～3个波传播的计算时间。伴随波场是通过反向传播伴随源,伴随源由目标函数关于合成数据的一阶导数给出。目标函数的一阶扰动为:

(2)

由此可确定伴随源s(t)=W(t)[p(t)-d(t)]。确定了伴随源之后,即可通过正传波场与反传波场的相关求取模型更新的梯度,继而通过不同的迭代求解算法逐步更新模型。

1.2 最优化传输函数的模式匹配特性

最优传输函数基于模式的匹配。这个性质可以用插值的例子来说明。图1在左下角显示圆盘,在右上角显示正方形盘。传统的插值(线性插值)将它们与两个加权因子分别相乘,之后相加。图2中的第1行显示了线性插值的演化,其中圆盘上的加权因子从1逐渐减小到0,同时正方形盘上的加权因子从0增加到1。可以看到,插值的效果是此消彼长的变化。但有时,希望可以像水流一样把物体形态光滑地变化成为另一个。图2中的第2行显示的正是这种变化,这一种基于最优化传输的插值,首先求出圆盘与正方形盘直接的传输平面,然后对平面插值即可。可以看到,圆盘的形状和位置均平滑地演化成正方形盘和其所在的位置。

最优传输函数的这种插值特性对FWI非常有帮助。传统的FWI经常被周跳问题困扰。最小二乘函数只是点对点地比较两个地震道。一旦两个地震事件相距较远,在伴随场震源上两个事件就不重叠,也就是说,它们之间就没有相互作用了。这就解释了为什么最小二乘函数只在两个地震事件误差在半周期内时起作用,同时在误差大于半周期的时候会陷入局部极值。而波动方程层析建模的凸性较FWI更好,因为两个事件总是可以通过互相关联系起来。因此,我们期望具有全局对比功能的最优化传输函数,它比最小二乘函数具有更好的凸性。

1.2.1W2函数

(3)

应遵守质量守恒原则,即:

(4)

图3 计算W2的示意

基于方程(3)和方程(4),可以计算出伴随源。伴随源s(t)由W2相对于合成的p(t)的偏导数给出[47]。引入中间变量,有:

(5)

最后的伴随源为:

(6)

数值实现过程可总结为:

1) 对每个地震道进行预处理,使其为正。建议将每个处理过的地震道除以其自身的质量,使总质量为1。

移动和归一化属于线性变换,把数据变为正信号还包括非线性变换。非线性变换包括平方、指数或双曲正切(tanh)。CHEN等[47]在地震定位的研究中采用平方法对地震数据进行预处理。作者设计了一个两层模型,并根据模型层位速度的不同计算了相对应的合成地震图,最终做了目标函数的凸性分析,如图4 所示。令人惊讶的是,将地震数据进行平方这个预处理优于移动和归一化策略。采用平方预处理得到的目标函数分布图更加光滑,凸性更好。

图4 目标函数误差分布[47]

与此同时,体波的振幅和弱事件的振幅会得到增强。振幅平衡是弹性全波形反演中重要的步骤之一,对于提升体波对于深部的反演贡献巨大。HE等[31-32]在各向异性介质中详细论述了振幅问题。

1.2.2BFM函数

一维问题的最优化传输存在解析解。与一维问题不同,高维问题不存在最优化传输的解析解。多年来,数学家们一直在努力寻找有效的求解方法。来回跳跃法BFM属于有效求解的方法之一。首先,对原问题构造拉格朗日函数,将其转化为对偶空间的等价问题。至此,有两种等价的最优化传输问题,一种在原始空间,另一种在对偶空间。BFM的独特优点是,最优解是在原始空间和对偶空间之间迭代产生的。在原始空间迭代一次,然后跳到对偶空间,在对偶空间同样只迭代一次,然后再跳回原始空间;一直循环直到问题收敛。这种方法的数学基础比较复杂[46]。BFM仍然要求输入是非负的,用于W2的预处理也适用于BFM。

1.2.3KR-OT函数

KR-OT是基于1-Wasserstein的函数;但并不是求解原始问题,而是转化到对偶空间中求解。KR-OT的代价函数是基于1范数c(x,y)=|x-y|,在对偶空间中的表达式为[30]:

(7)

式中:φ(x)在有界的1-Lipschitz空间中,而λ是一个用户定义的参数。第1个不等式保证了φ不能太大的要求,第2个不等式排除了突然变化的情况。HE等[32]仔细研究了不等式约束的影响,其结论是存在较大的λ取值空间,都可以给出合理的解决方案φ。KR-OT的显著特性是,该公式只依赖于两个分布之间的差异,所以不要求输入的地震信号是正的。对于任何负分布或非正分布,可以将两个分布加上同一个常数变为正值,但是从KR-OT的表达式看,这个常数是会被消去而不起任何作用。

在全波形反演中,地震数据是以离散形式存储的。在三维情况下,对应的KR-OT公式为:

(8)

式中:常数因子Δx1,Δx2,Δx3在公式(8)中的第1个方程中被省略,uijk是合成数据与观测数据的差异。公式(8)可以用SDMM(simultaneous direction method of multipliers)高效求解。METIVIER等[30]详细说明了如何解决KR-OT。其中,求解矩阵逆是一个沉重的计算负担。后来,METIVIER等[48]证明了其中的矩阵运算与泊松方程是等价的,从而可以在频率域用FFT方法或在时域用多重网格方法有效地求解泊松方程。虽然如此,三维计算仍然非常昂贵,主要原因是每次求解KR-OT都需要多次求解三维泊松方程。在二维空间中,从公式(8)中去掉一维,可以大大减少计算时间。因此,伪三维处理也是折衷的方法。将三维地震数据体沿地震数据线方向分割成切片(地震数据切片形成二维矩阵),并使用2D KR-OT求解器对每个切片进行处理。

KR-OT的输入信号是残差,故所有的地震事件从一开始就是混合的。这会影响KR-OT的凸性。与其它OT公式相比,KR-OT的独特优势是KR-OT将解约束为1-Lipschitz函数。这种优势对于包含面波的陆地地震数据具有很大的潜力。图5为KR-OT能平衡振幅,展示的是一个陆地数据全波形反演的伴随场。L2伴随源以面波为主,因此,重建由体波采样的深层区域是困难的。相比之下,KR-OT伴随源的振幅分布平衡,反演确实显示了从浅到深的重构模型。

图5 KR-OT能平衡振幅[32]

1.2.4GS-1D函数

图空间最优化传输(Graph Space Optimal Transport)处理的是离散化的问题。将一维地震道扩展为时间-振幅的二维空间,最优化传输的代价函数为:

cij=(pi-dj)2+η(ti-tj)2

(9)

式中:pi是在ti时刻采样的合成地震道;dj是在tj时刻采样的观测地震道;cij是点(ti,pi)至点(tj,dj)的移动成本;η是用户定义的参数,用来控制时间轴的权重。最优化传输要求质量守恒,在图空间中,质量守恒约束指的是合成地震道与观测地震道包含相同数量的点。配对一对地震道的总成本是对所有个体成本的求和,即:

(10)

(11)

伴随源si看起来像传统的L2的情况,但其中的观测数据是重排列之后的观测数据。在L2的情况下,两个不重叠地震事件是没有相互作用的,这也是周跳的根源。GS 1D依然能够通过求解最优的排列来联系这些事件。

图6显示了L2、KR-OT和GS-1D的伴随源。图6a中,观测数据包含多个地震事件,而合成数据只有两个。图6b中,常规的最小二乘函数无法平衡振幅,而KR-OT和GS-OT都在很好的平衡振幅。输入信号由Ricker子波组成,模拟的是事件错位和事件缺失的情况。L2伴随源具有较大的振幅差异,而OT伴随源确实平衡了振幅。然而,这并不意味着在这两个OT函数中都忽略了振幅信息。在全波形反演中,随着迭代的继续,观测数据的振幅会得到越来越好的利用。最后,所有的事件都是完全拟合的。

图6 L2、KR-OT和GS-OT的伴随源[31]

1.2.5GS-2D(高维OT)

(12)

图7和图8均引自XU 等[45]。图7显示了两个二维质量分布之间的映射;两个质量分布都是由高斯函数生成,但是中心位置不同。二维的GS-OT算法很好地求解出了匹配关系(黑色实线)。若是采用一维算法,则只能沿着横向或者纵向匹配,而无法正确求解出真实的匹配关系。在全波形反演中,地震事件经常是在时间轴方向和空间轴方向都有移动,这也是为何要发展高维OT算法的原因。图8显示的是三维质量分布之间的三维映射。两个球体分别对应的是两个三维的质量分布,可以看出,三维的GS-OT算法正确地匹配到了两个三维质量分布。XU 等[45]将高维最优化传输GS 2D应用于声学和弹性全波形反演试验,包括弹性SEG/EAGE逆冲模型和Chev-ron 2014的数据盲测试验。在逆冲模型中,L2反演重构出了横波速度参数vS的浅层结构,GS 1D重建了vS的浅部和深部。相应地,体波和面波都比L2得到了更好的利用。GS-2D得到了最佳的反演,即vP和vS都显示了高分辨率和高保真度的结构。在XU等[45]展示的Chevron 2014数据盲测实验中,使用熵变换的代价函数x=xlogx得到了更干净的速度。

图7 二维质量分布之间的最优化映射

图8 三维质量分布之间的最优化映射

2 反演实验

我们已经实现了上述所有的最优化传输函数W2,BFM,KR-OT,GS-1D和GS-2D的FWI。本节将展示两个模型测试结果。第1个是透射例子;第2个是选取的SEAM Ⅱ Foothill模型的切片。

2.1 透射模型

真实模型是在均匀背景vP=3.0km/s的模型中添加了vP=3.6km/s的圆形异常体。圆形异常体的直径为3.6km。假设密度为常数ρ=1g/cm3。以均匀背景速度作为初始模型。图9显示了真实的模型和初始的模型。在顶部等距离布置了32个爆炸源,在底部布置了300个压力接收器。震源子波为10Hz的Ricker子波。以10Hz为参考计算,圆形异常体的大小是主波长的10倍。

图9 全波形反演测试

图10显示了观测数据以及观测数据与合成数据的交错比较。对比图清楚地说明了周期跳跃问题。L2反演进行了6次迭代陷入了局部极值,反演结果在图11a中。KR-OT反演也被局限在某个局部最小值中(图11d)。在这个特定的例子中,初始模型与实际模型相差太远。正如作者所解释的,KR-OT将L2误差作为输入,它的性能是受到影响的。其它OT反演(W2,BFM,GS-1D,GS-2D)都正确地恢复了模型。

图10 观测数据、合成数据与观测数据对比

2.2 SEAM Ⅱ Foothill模型

所用模型是从原始的3D SEAM Ⅱ FOOTHILL模型中抽取的一个剖面,同时去除了近地表低速区。真实模型如图12所示,纵、横向采样均是20m,横向有720个格点,长约14km;纵向有200个格点,深约4km。可以明显看到地层发生扭曲和位移,但是模型沿着纵向的整体变化率并不大。采用3Hz的Ricker子波作为震源子波。在顶部布置了48个震源和270个接收器。模型正演采用了吸收边界条件。图13为观测数据的炮集。图中的初至大致沿着一条直线,因为模型纵向变化不大,所以潜水波(diving wave)穿透深度较浅,从而导致走时跟距离接近线性关系。本文创建了4个初始模型,如图14 所示。前3个由真实模型采用高斯平滑得到,平滑半径分别为10,20,30个网格长度。最后一个是一维线性递增模型。相应的合成地震数据如图15所示。初始模型越简单,合成数据也越简单。对前面两个初始模型,尚能识别反射波;后面两个不易辨认反射波。对每个初始模型都进行了6次反演测试,包括1次常规反演和5次OT反演。总共用SEAM Ⅱ Foothill模型进行了24次反演实验。

图12 实际速度模型

图13 观测数据

图14 初始模型

图15 根据图14计算出来的合成地震数据

第1个初始模型最接近真实模型,且其初始合成数据已经包含了大量反射事件。图16显示了使用第1个初始模型的反演结果。所有6个反演都给出了正确的结构,且反演结果基本相当。常规全波形反演中的L2目标函数的吸引域(attraction basin)较小,从另一个角度来说,L2的分辨率较高。图16的反演结果显示了OT函数也同样具有较高的分辨率。与此相比,若是采用波动方程层析(wave equation tomography),通常需要用L2接着往下反演;而OT则具有一步到位的潜力。

第2个初始模型比第1个更为光滑,大部分清晰的反射结构已经被平滑掉。其合成数据仅能识别一个大的反射事件。初至波则与第1个模型的初至波较为接近。原因为:①初至波穿透较浅;②初至波沿着浅地表传播,其到时本身就是模型结构的加权平均,所以对模型光滑不会导致初至太大的变化。图17 展示的是6个目标函数对应的反演结果。常规反演没有完全重构速度模型,其最深处能识别的结构在中间2km深处。与真实模型相比,其浅部的速度缺乏背景场的更新,更接近于偏移模式。在早期的全波形反演发展中,地球物理学家已经意识到反射波对于中间波长的速度模型是缺乏反应的(intermediate-wavelength gap),所以用L2函数会导致全波形反演向着偏移的模式演化,并陷入局部极值。图17所示结果实际上揭示了OT函数对于走时信息的提取能力。OT函数是基于模式匹配的,模式上的匹配是隐含着运动学信息(kinematic),如图17和18所示。一维OT提取的是走时信息,高维OT除了提取走时信息外,还提取了事件沿着空间轴方向的移动信息。这些运动学信息的提取对于模型低波数成分的恢复起了决定性作用。从图17的5个OT结果上看,3个二维OT函数的反演结果略微优于两个一维OT函数的结果。在模型的左下方,两个一维OT的结果未能清晰地重构出低速带。需要说明的是,一维GS-OT的用户输入参数不大容易选取最优输入参数,在特定模型上累计的经验推广到另一个完全不同的模型上往往并不会得到同样的结果。从原理上讲,GS-OT是具有良好的凸性的,比W2更优越,图17e中展示的反演结果或许还有进一步提高的可能。

在作者积累的全波形反演经验中,模型陷入局部极值后可能会有一些难以解释的现象。图18展示的是用第3个初始模型得到的6个反演结果。相比于图17中的常规反演,图18中的常规反演得到的模型更新更深。虽然两者都是因全波形反演的偏移模式而陷入局部极值,但无法解释为何图18a会有更深的模型更新。与图17类似,图18中同样是3个二维OT函数的反演结果略微优于两个一维OT函数的结果。并且GS-1D的反演结果确实比W2更好。

图18 以图14c作为初始模型得到的反演结果

图19展示的是初始模型为速度随深度线性增加模型对应的6个反演结果,对于所有这些函数,线性增加的模型可能离真实的模型太远。W2的反演可能是因为在陷入局部极值的情形下步长不合适而出现了特别大和特别小的速度结构。其它5个反演结果应当是全波形反演的偏移模式,背景模型难以更新。值得注意的是,KR-OT对应的反演结果模型迭代得最深,这可能跟KR-OT独特的平衡能力相关。KR-OT的输入信号即是数据残差,在经过KR-OT平衡后,来自于深部的反射信号得到了加强,所以导致较深的模型迭代。

图19 以图14d作为初始模型得到的反演结果

3 讨论

图11中的常规全波形反演陷入局部极值,即使没有真实模型参考也能识别是局部极值。但在实际的全波形反演中,应用常规全波形反演后遇到的情形并不像图11那么清晰。关于实际数据中出现的局部极值现象,一个更为准确的推断是,在部分区域模型恢复得很好,部分区域与真正的结构差异较大。所以导致了全波形反演与某些测井资料对应得很好,而与其它测井信息则差异较大。也正是因为这个原因,才凸显出凸性更好的目标函数的重要性,这样即可在绝大部分区域都能成功避免局部极值。从目标函数的研究历史看,从最初的点对点的比较,逐步发展为道与道的比较。基于道比较的目标函数包括:包络函数、互相关、反卷积和互相关走时差等一系列函数。最优化传输是基于道的整个数据体的模式比较,支持整个数据体之间的比较。从研究历史看,基于数据整体的比较是一个发展方向,其优势在于可以在全局比较中正确匹配地震事件的偏差,是全局优化策略。

最优化传输函数的凸性取决于:①范数。通常会选择1或2范数(p-Wasserstein函数)。通常可以认为2范数比1范数更凸。此外,在大多数情况下,2范数更容易操作。相比之下,1范数在零点附近不是平滑的,比如,|x|在零附近的导数就是间断的。不恰当地将1范数合并到全波形反演中可能会导致伴随源在某些点或者区域的不连续性。KR-OT函数虽然是基于1范数,但并没有出现这种不连续的现象,因为KR-OT是在对偶空间中表示的,1范数被转移到了约束方程中。②预处理。对于数据预处理需要谨慎对待。比如W2常用到的移动加规则化,其移动量不宜太大。基于非线性变换的预处理有很多,在作者的经验中,tanh表现不错,而且用户输入参数容易调节。③时间空间坐标是否引入最优化传输目标函数,图空间最优化传输是将时间坐标和空间坐标显式地引入,其结果也变得更凸,同时,避免了负信号的处理。

具有全局优化潜力的最优化传输函数的计算量比常规函数要大。在实用化方面,一维的最优化传输函数的计算量很小,可以应用在二维和三维全波形反演中。若将三维数据分割成二维的矩阵,那么,本文提到的3个二维最优化传输函数同样可以应用,计算代价在可接受范围内。若是直接处理三维数据体,计算量较大,可能需要一些加速策略。从目标函数的性能来看,没有哪一个函数可以处理所有的情形,但是有些函数会在大部分地质环境中表现良好。使用SEAM Ⅱ Foothill模型的全波形反演实验证明,总体趋势是2D最优化传输函数优于1D版本。BFM功能优于W2功能;GS-2D技术优于GS-1D技术。从理论上说,二维函数能提取空间轴方向的信息,比一维更具潜力。这与作者在其它各种模型上测试得到的认识是一致的。综合作者在各种模型上的测试经验看,W2在数据模式较为简单的情形下表现良好,可能是因为预处理影响到了函数性能。KR-OT因为输入数据是数据误差,导致地震事件从一开始就混杂在一起,所以合成数据跟真实数据不能相差太远。但是其具有独特且优越的振幅平衡能力。在陆地数据的弹性全波形反演中KR-OT是首选。从原理上说,KR-OT比L2更凸,但会比GS-1D略逊一筹。GS-OT(包括GS-1D和GS-2D)的优势在于将坐标引入计算代价中,从根本上避免了数据的预处理。但是,凸性良好的函数具有共同的小瑕疵,即在反演模型接近真实模型时,反映在数据上的偏差可能不会太大,所以反演收敛速度会很慢。基于此,在GS-OT反演后,可以酌情在后续加入常规反演。BFM的性能原则上比KR-OT更凸,但逊于GS-2D。BFM是基于2范数的最优化传输,而KR-OT是基于1范数的,2范数比1范数更凸。GS-2D也是基于2范数,同时,避免了数据预处理,应当比BFM更凸。

4 结论

本文综合对比了目前学术界和工业界新提出的最优化传输目标函数。常规的最小二乘目标函数是点对点的比较方式,对于勘探地震数据地震道中的时间相关性和地震道之间的空间相关性没有显式地提取。当初始模型足够准确,合成地震事件与真实地震事件偏差在半个波长以内的时候,常规的最小二乘目标函数是可以反演出足够接近真实的地下速度模型。然而在实际数据处理中,相对于采集中设定的低频范围,初始模型常常是不够准确的。这也是学术界和工业界寻求其它目标函数的动力。从目标函数的研究历史看,的确是朝着高维方向发展。在目前的研究中,最优化传输通常采用1范数和2范数;对于给定的范数,又可以分为单道处理和高维处理。综合我们的测试结果看,高维的最优化传输函数优于一维的版本。在实际应用中,可以综合考虑计算成本和效果。对于工区不太复杂的情形,可以采用W2,因为其存在解析解,计算速度最快。在复杂地区,可以采用BFM,KR-OT,GS-1D或者GS-2D;在这4种最优化传输函数中,以GS-1D计算成本最小。其中,只有KR-OT是基于1范数,其凸性可能逊于BFM;但是KR-OT独有的振幅平衡能力使得KR-OT在处理存在强振幅差异地震数据的时候具有特殊的优势。