四元数在多分量地震数据处理中的应用研究现状

李远芳,汪 超,王 赟

(1.中国科学院地球化学研究所,矿床地球化学国家重点实验室,贵州贵阳 550081;2.中国科学院大学,北京100049;3.中国地质大学(北京)地球物理与信息技术学院,北京 100083)

多分量地震勘探技术利用了地震波的纵横波特性。检波器记录到的多分量地震数据保存了质点在垂直和水平方向的线性运动情况,完整记录了地震波的矢量场信息[1]。矢量场信息有助于更加详细、精确地刻画地下介质的构造、岩性、流体饱和度、孔隙压力和裂缝等特征[2]。目前,处理多分量地震数据通常采用先将矢量场信息分解为几个标量场分量数据,再对各个标量场分量数据单独处理的方法[3]。基于标量场的多分量地震数据处理技术虽然可以处理多个分量数据,但因未能充分挖掘多分量地震数据携带的矢量信息,容易破坏矢量数据的相对振幅关系[4-8]。因此,保持矢量特征的矢量场处理技术是当前多分量地震数据处理技术发展的重要方向之一[9-12]。

1840年,哈密顿首次提出四元数概念,上世纪60年代,四元数广泛应用于刚体力学等领域,而后经过几十年的发展形成了较完备的四元数理论。理论上,四元数作为Clifford代数的一个子代数,是复数在三维实空间的延伸。从形式上看,四元数由一个实部和3个虚部组成,适合表达多分量数据,其在描述和分析处理三维空间各种矢量运动上具有独特的优势[9]。四元数的建立为矢量信号处理提供了有力的工具,将四元数理论推广至多分量地震数据处理领域,可以建立一个完整的矢量场多分量地震数据处理流程[13-17]。1989年,DAVIES等[10]首次系统提出了地球物理积分变换中的四元数方法,该方法在二维地球物理场方程计算结果的基础上引入四元数,为三维向量场的积分变换提供了有效的手段。将四元数与数学变换工具相结合,可以使得一些标量场的处理手段能够被应用于矢量场的处理中,如CHOUKROUN等[11]和SABATINI[12]将卡尔曼滤波器(Kalman filter)推广到四元数领域并将其应用于磁传感器方向定位。

1 四元数基础理论

1.1 四元数及四元数矩阵

四元数通常采用下列形式来表示[9]:

q=a+bi+cj+dk

(1)

式中:a,b,c,d均为实数,i,j,k为虚数单位。i,j,k满足如下关系:

ij=-ji=k

jk=-kj=i

(2)

ki=-ik=j

当a=0时,(1)式为q=bi+cj+dk,q为纯四元数。

四元数是复数的扩充形式,满足一般复数的运算法则,但四元数乘法不满足乘法交换律。为了避免乘法的不可交换性造成的繁琐运算,通常引入四元数复表示和实表示将四元数乘法转化为复数域和实数域上的等价运算问题。定义四元数q=a+bi+cj+dk的实表示为:

(3)

四元数q的复表示为:

(4)

同样地,定义四元数矩阵A如下:

A=(aij)m×naij∈

(5)

(6)

(7)

1.2 多分量地震数据的四元数形式

四元数形式上的特点使得它非常适合表达多个分量的数据。例如海底电缆采集的四分量(4C-OBC)地震数据,即地震检波器记录到的压力标量和x,y,z3个方向上的位移矢量,可以表示为四元数形式,3个方向上的位移矢量对应3个虚部,压力分量对应实部。而陆地勘探的三分量地震数据,只需将实部设置为0,表示为只含有3个虚部的纯四元数即可。

若用pi(t)代表第i时刻记录的压力分量时间序列,用xi(t),yi(t),zi(t)分别表示第i时刻地震检波器记录的X,Y,Z分量上的时间序列,则可以定义四分量地震数据时间序列如下:

qi(t)=pi(t)+ixi(t)+jyi(t)+kzi(t)

(8)

类似地,对于三分量地震数据来说,可以将其定义为一个纯四元数的时间序列并表示如下:

qi(t)=ixi(t)+jyi(t)+kzi(t)

(9)

如果频率域两个分量地震数据用D1,D2表示为:

D1(ω,x)=D1,real(ω,x)+D1,imag(ω,x)i

(10)

D2(ω,x)=D2,real(ω,x)+D2,imag(ω,x)i

(11)

那么二分量地震数据可以用四元数形式表示为:

Q(ω,x)=D1,real(ω,x)+D1,imag(ω,x)i+D2,real(ω,x)j+D2,imag(ω,x)k

(12)

一个大小为m×n(m为时间,n为道数)的多分量地震记录可以统一表示成m×n大小的四元数矩阵:

(13)

因此多分量地震数据可以表示为四元数矩阵,并利用四元数的代数方法来处理。

四元数应用理论和方法的不断完善为利用四元数进行多分量地震数据矢量去噪研究提供了坚实的基础和新工具。MANDIC等[18]提出利用HR梯度算子对线性四元数进行求导运算。JIANG[19]对HR梯度算子进行扩展,提出了适用于非线性四元数的广义梯度算子,并给出了最小二乘法和非线性自适应法计算公式的导数表达式。XU等[20]推导了四元数优化问题中梯度、海森矩阵(Hessian)和学习进化的运算法则。由于四元数在矢量代数和矢量分析中的优越性,因此在信号处理领域中,相关学者将四元数理论与现有的信号处理手段相结合,发展了一套较为完备的高维矢量信号处理方法。1992年,ELL[21]首次给出了四元数傅里叶变换(quaternion Fourier transform,QFT)的定义,随后SANGWINE等[22]和ELL等[23]对离散四元数傅里叶变换进行了深入研究。2001年,PEI等[24]利用快速傅里叶变换算法(fast Fourier transform,FFT)实现了四元数傅里叶变换、四元数卷积和四元数相关计算。BAHRI等[25]提出了短时四元数傅里叶变换,可以对非平稳信号进行更加准确的分析。BÜLOW[26]在四元数傅里叶变换的基础上,提出了Gabor小波变换,实现了具有局部高斯窗的四元数傅里叶变换。在此基础上,SOULARD等[27]提出了可对四元数时变系统进行多尺度多分辨率分析的四元数小波变换。

2 四元数在多分量地震数据处理中的应用

2.1 速度分析

速度分析是地震数据处理流程中必不可少的环节之一,对后续的叠加成像至关重要。速度分析方法主要包括两类。①叠加类方法:对各个地震道的振幅值求和并计算平均振幅,最大振幅值所对应的速度即为叠加速度值;②相关类方法:计算各个地震道的相关系数,根据相关系数值来获取速度。叠加类方法计算量较小,计算结果可靠性高;相关类方法的灵敏度高,求取的速度谱峰值明显,但随机干扰较大时,其计算结果不稳定。

对于多分量地震数据,通常将其分为3个方向的单分量标量地震数据分别进行处理。单分量的标量速度分析方法计算过程简单但却忽略了如纵横波速度比、地震事件相关性等波场矢量信息。为了保留波场矢量信息,GRANDI等[14]引入四元数来实现三维数据的速度分析,即先将时窗数据映射到四元数域中,再将常规的速度分析方法应用于四元数域。四元数速度分析流程如图1所示。

图1 四元数速度分析流程

图1中的速度检测因子即为各种速度分析方法(即相似系数法、协方差法、匹配滤波法以及复匹配滤波与协方差结合法)所用的相干函数。上述速度分析方法及相关计算公式分别介绍如下。

2.1.1 四元数相似系数法[14]

该方法的相似系数Qs为:

(14)

式中:T是时间窗口的长度;G是地震检波器的数量;qi为(8)式中四元数形式的多分量地震时间序列。四元数相似性系数法需要在CDP道集上对双曲线时窗内的多道炮集记录进行互相关,互相关结果反映了地震信号能量的变化。相似系数作为最常用、最简单的相干函数,为所有地震道叠加后能量和叠加前地震道能量和的比值,其计算过程简单,但对于相近或干涉的同相轴,因叠加的能量存在极大值,会影响速度谱的分辨率。此外,该方法计算高效,鲁棒性强,但当振幅值变化较大时,计算得到的相似系数值不准确[28-30]。

2.1.2 四元数匹配滤波法[14,31]

利用该方法模拟的波形与给定波形的相关度为:

(15)

式中:τi为双程旅行时;vi为叠加速度;ρ(xj;τi,vi)为给定的单道归一化相关函数,ρ(xj;τi,vi)=w*d(xj)/‖w‖‖d(xj)‖,其中,w为已知的子波,d(xj)为经过NMO校正过的时窗数据;*代表褶积;‖·‖代表范数。

匹配滤波法建立在已知地震子波的振幅谱基础上,根据不同速度的分量频谱不同的原理,对每个分量使用不同的地震子波。单分量数据处理时,对时窗内的数据进行希尔伯特变换得到复数域的矩阵,而四元数矢量数据处理时直接将三分量地震数据映射到四元数域得到四元数矩阵。根据(15)式对四元数矩阵进行匹配滤波,能有效压制随机噪声和相干事件的影响。与四元数相似系数法相比,四元数匹配滤波法可用于处理振幅随偏移距变化的地震资料。

2.1.3 协方差方法[14]

该方法主特征值与其它特征值的信噪比SNR为:

(16)

式中:λi为协方差矩阵的特征值;N为总道数。

GRANDI等[14]将协方差方法应用于超复数域四元数中,首先将每个地震道表示成四元数形式,给定一组速度范围值,从各道对应的速度时窗中提取四元数矩阵,然后,计算提取的四元数矩阵的转置矩阵,四元数矩阵与四元数矩阵的转置矩阵相乘以构建一个四元数协方差矩阵,再计算协方差矩阵的特征值,最后利用(16)式计算出的信噪比值来选择权重得到速度谱。当拾取速度正确时,信号的能量只反映在主特征值中,其它的小特征值急剧减小。对多道求和得到的超道集采用上述方法计算可以提高信噪比,减少计算时间。

2.1.4 四元数匹配滤波和协方差结合法

该方法将四元数匹配滤波和协方差方法相结合,进一步提高了时间分辨率,同时有助于消除不同接收器响应对地震数据的影响[15]。

以上4种方法在合成地震数据以及实际地震数据上的测试应用均表明四元数相干函数拾取到的速度谱比对应的标量相干函数拾取到的速度谱具有更高的分辨率,同时拾取时间比X,Y,Z每个分量逐一拾取的总时间更少,拾取速度更快[14-15,32-33]。更为明显的优势在于矢量数据的保真度高,对于存在明显的剪切波分裂的情况,四元数相干函数可以同时聚焦快横波(S1)和慢横波(S2)地震事件。

2.2 反褶积处理

反褶积处理是地震数据处理的基本步骤之一,其主要目的是压缩地震子波,提高地震资料纵向分辨率。多数反褶积方法基于标量场,对矢量场多分量反褶积的研究较少。1970年,TREITEL[34]设计了多通道反褶积,随后CLAERBOUT[35]将其发展为多分量反褶积,他们将所有分量的时间序列排列为一个长矢量,作为一个单分量地震记录进行处理。2012年,MAZZOTTI等[15]和MENANNO等[36]利用四元数表示矢量的优势,提出了反褶积的四元数形式,并且推导了维纳反褶积的四元数公式,将经典的维纳反褶积方法拓展为四元数维纳反褶积方法。合成地震记录和实际地震记录中标量反褶积和四元数反褶积处理结果表明,四元数维纳反褶积方法无论应用于脉冲反褶积处理还是预测反褶积处理,都优于标量反褶积方法。

四元数维纳反褶积方法是对复维纳反褶积方法的拓展,设计四元数滤波器fn使得实际输出yn和期望输出dn之间的平方差误差最小,其中实际输出取决于对输入道xn进行的四元数滤波效果。误差函数J为:

J=E[|dn-fn*xn|2]

(17)

式中:dn为期望输出;fn*xn为实际的输出;E为期望。

MENANNO等[36]将复数求导算子推广到四元数中,利用对误差函数求导来求取误差函数最小值。由于四元数乘法具有不可交换性,所以四元数维纳滤波的矩阵形式为:

(18)

式中:Rxx[i],Rdx[i]分别为四元数自相关和互相关函数,分别定义如下:

(19)

(20)

选取相同的参数分别进行四元数矢量反褶积和标量反褶积处理,处理结果表明水平分量四元数矢量反褶积处理方法明显优于标量反褶积处理方法,其输出更接近理想反射系数序列。垂直分量四元数矢量反褶积和标量反褶积处理结果差异不明显,但根据每个分量实际输出和期望输出之间的平方误差值可知,多分量四元数矢量反褶积处理结果更接近期望输出。为测试不同反褶积方法的处理效果,对同一个合成地震记录加上80m深水层引起的海底混响,分别进行X分量四元数矢量反褶积和标量反褶积处理,图2,图3和图4分别为X、Y、Z分量四元数矢量反褶积处理和标量反褶积的处理结果及二者的误差平方和,图中的红线代表标量反褶积的处理结果,绿线代表多分量四元数矢量反褶积的处理结果,可以看出,对水平分量(X和Y)而言,四元数矢量反褶积处理结果明显优于标量反褶积处理结果。因Z分量四元数矢量反褶积和标量反褶积的信号几乎完全重叠,故对二者的误差平方和进行了计算,计算结果表明Z分量四元数矢量反褶积处理的误差平方和小于标量反褶积处理误差平方和。

2.3 矢量去噪与矢量重建

基于矢量场的多分量地震去噪方法[37-48]包括3种方法。①极化率滤波方法:根据不同类型的地震波在传播过程中的偏振方式不同进行差异识别和波形分离。②基于矢量统计排序的矢量波场滤波方法:作为一种以多元数据排序统计分析为基础的滤波方法,可细分为矢量中值滤波方法、矢量方向滤波方法和多道α截集均值滤波方法。③基于组稀疏表示的多分量联合去噪方法:将时间域的含噪信号通过数学变换到其它域进行信号与噪声的分离,其中的数学变换通常包括傅里叶变换、小波变换和S变换等。利用张量降秩理论中经典的奇异值分解方法可以实现四元数多分量地震数据去噪。

图2 X分量四元数矢量反褶积处理和标量反褶积处理结果(a)及二者的误差平方和(b)

图3 Y分量四元数矢量反褶积处理和标量反褶积处理结果(a)及二者的误差平方和(b)

图4 Z分量四元数矢量反褶积处理和标量反褶积处理结果(a)及二者的误差平方和(b)

2004年,BIHAN等[13]首次将四元数应用于多分量地震数据处理中,提出了一种基于四元数奇异值分解的去噪方法。该方法首先将3个分量的地震数据分别作为四元数的3个虚部,而后将多道地震数据构成一个四元数矩阵,再对四元数矩阵进行奇异值分解,最后选取能量集中的前几个奇异值进行重构达到去噪的目的。四元数奇异值分解算法如下:

1) 给定四元数矩阵A,求出A的复表示矩阵Aσ;

2) 对复矩阵Aσ进行双对角化分解,得到双对角化矩阵B;

3) 对B进行奇异值分解,求出Uσ和Vσ;

5) 利用复表示矩阵的逆变换,求出U,V,Σ。

图5,图6,图7和图8分别是原始的不含噪三分量地震信号、含噪的矢量地震信号、四元数奇异值分解去噪后的矢量地震信号以及标量奇异值分解去噪后的矢量地震信号,3个分量分别用四元数的3个虚部单元i,j,k表示。不难发现,四元数奇异值去噪方法比标量奇异值分解去噪方法能更好地压制噪声且保幅效果良好。

去噪理论的核心思想为放大信号和噪声之间的差异,通过识别差异来保留需要的信号。基于张量降秩理论的去噪算法原理为,一个张量或低秩的矩阵所表示的完整的地震数据,会因噪声增加矩阵的秩或张量的阶,故利用降秩可以有效地对地震数据进行去噪。缺失的地震道被视作随机噪声,影响矩阵的秩,因此降秩算法能够同时实现去噪与重建。基于稀疏变换理论的去噪方法如QFT、Radon变换和Curvelet变换[37]方法等均可用于地震矢量数据的重建。STANTON等[17]将四元数与稀疏变换理论相结合,提出了一种将四元数傅里叶变换与凸集投影(projection onto convex sets,POCS)相结合的多分量地震数据矢量重建方法,首先将X和Y两个水平分量变换到频率域以获得4个正交量(每个分量均含1个实部和1个虚部),然后将4个正交量用一个形如(12)式的四元数统一表示,对每个频率成分沿4个空间方向进行QFT,再将正交的X和Y两个水平分量数据变换到频率波数域,最后利用POCS方法重建多分量地震数据。该方法保证了输入分量的正交特性,同时提高了重建地震数据的质量,但该方法存在局限性,其在多分量数据频谱不重叠时重建地震数据的质量较差。

图5 原始的不含噪三分量地震信号

图6 含噪的矢量地震信号

图7 四元数奇异值分解去噪后的矢量地震信号

图8 标量奇异值分解去噪后的矢量地震信号

2.4 其它应用

白清云[49]首先利用二阶对称的四元数QFT在一个象限内完全恢复了原信号,然后定义了二维的四元数解析信号,再对四元数解析信号反演求取纵横波的二维瞬时属性比值,得到的油气藏刻画结果更清晰,最终能更有效地进行油气预测。从四元数信号中提取的瞬时振幅比属性充分考虑了储层横向变化带来的影响,因此油气预测效果优于利用标量信号提取出的纵横波瞬时振幅比属性。四元数还可应用于地球物理勘探中的采集环节,KRIEGER等[50]提出了一种基于四元数的多分量地球物理传感器优化定位的方法,利用该方法计算得到的多分量检波器方位矫正量更加准确稳健。SAJEVA等[16]提出了一种基于四元数奇异值分解(quaternion singular value decomposition,QSVD)识别体波与面波并从波场中分离瑞利波的方法。单色瑞利波模态的位移可以表达成秩为1的四元数矩阵,通过群速度校正、瞬时极化校正以及逐道相移旋转,可以使得修正后的瑞利波模式映射到四元数体中的第一个特征值中,因此采用QSVD方法能从体波和高阶瑞利波模态中分离出瑞利波。合成数据测试结果表明,采用QSVD方法可以分离体波与瑞利波,可以从高阶瑞利波中将基阶瑞利波分离出来。

四元数还可以应用于其它多分量地震矢量场中,如时间延迟分析和边缘检测[51]以及各向异性反演等[52]。

3 讨论与结论

四元数理论并非单一的理论,而是一个理论框架,将四元数方法理论与多分量地震数据处理相结合,以四元数矩阵傅里叶变换、四元数卷积、四元数矩阵奇异值分解等理论作为主要数学工具,辅以其它信号处理方法,如主成分分析、维纳滤波、采样定理等,能保幅保真地对多分量地震数据进行矢量处理。四元数背景下的一系列数学变换和信号处理理论的不断发展,为地球物理领域的多分量数据的矢量处理提供了新的研究方法和研究方向。四元数理论在多分量地震数据矢量处理中的可行性主要包括以下3个方面。

1) 四元数的四分量形式非常适合描述多分量地震数据,OBC采集的地震数据可以用实部来描述压力分量,虚部来描述3个位移分量,而陆地采集的三分量地震数据也可以用纯四元数来表达。

2) 四元数的本质是超复数,适于处理矢量数据,可以更好地保持各分量之间的相对振幅关系不变,保证矢量波场不畸变。从现有的应用可以看出,它充分考虑了地震波场的矢量性质,因此可以同时处理多分量地震数据并且直接以矢量形式输入。

3) 四元数相关理论发展较为成熟,相应的数学变换工具也较为完备。因此可以将标量场已有的地震数据处理方法移植到四元数域中,在四元数域中进行多分量地震数据的矢量处理。前述的四元数速度分析、反褶积、矢量去噪及矢量重建均采用标量场数据处理方法,引入四元数可以将已有的标量处理方法拓展到四元数中,进而推导出适用于四元数的计算公式。

虽然四元数在多分量地震数据矢量处理中具有一定的效果,但尚未形成一套完整的处理流程,四元数只零散地分布于多分量地震数据处理流程的少数独立环节上,其前后配套的处理技术仍不完善,目前可以从以下3个方面展开深入思考和研究。

1) 建立基于四元数的多分量地震数据的完整矢量处理流程,包含去噪、速度分析、动静校正、叠加、偏移等。

2) 在目前四元数多分量地震数据矢量处理技术的基础上,利用四元数研发或改进现有的标量处理技术。基于四元数奇异值分解的去噪方法只能压制水平同相轴的噪声,还不能处理倾斜同相轴的噪声,可以考虑采用四元数Cadzow滤波、四元数奇异谱分析法等对矢量地震信号进行去噪处理。此外,可以考虑将四元数拓展至其它经典的去噪理论,如四元数Kalman滤波等。

3) 以上四元数的应用仅限于人工源地震勘探领域,未来应考虑四元数在天然地震中的应用,如震相拾取、地震监测及定位等。

综上所述,四元数理论应用于多分量地震数据处理是可行且有效的。基于四元数理论提出更加完善、系统的多分量地震数据矢量处理方法,可以充分挖掘多分量地震信号携带的信息,该方法具有解决复杂各向异性介质问题的潜力。