2. Modeling and Imaging Laboratory, University of California, Santa Cruz, U.S.A. 95064
2. Modeling and Imaging Laboratory, University of California, Santa Cruz, CA, U.S.A. 95064
地震偏移成像技术中,各种局域化传播算子在处理速度横向快速强变化介质时显示出明显优势.局域传播算子一般是指空间和传播方向上的局域化,空间的局域化使传播算子更容易适应介质在空间上的变化,传播方向的局域化使波通过介质后演化相对简单,从而能够保持较高精度.这些方法中,均采用一组局域化基本函数来表示波场,然后通过延拓这些基本函数来达到波场外推的目的.这些基本函数在介质中传播时必须遵循波的传播规律,可以按波的高频近似方法直接从地表传播到目标区(渐近方法),也可以按波动方程每次传播一个小的步长分多次来完成(波场步进延拓方法).本文主要讨论基于单程波动方程的局域化波场延拓方法.
Steinberg等[1-2]首先以高斯函数为窗函数,使用加窗傅里叶变换和扰动理论,推导了局部传输算子,并给出了波在局部角度域传播的步进算法.然而,加窗傅里叶变换分解和重建非常费时,限制了该方法在实际中的推广应用;另外,该方法对偏微分方程使用全局扰动求解,仍然不能适应横向快速强变化介质.为了充分发挥局域化传播算子的优势,Wu等[3]提出了基于局部背景参考速度和局部扰动理论的小波束(beamlet)偏移方法,其中主要发展了基于Gabor-Daubechies标架[4-5]和局部余弦变换[6]的小波束偏移.Gabor-Daubechies标架分解与加窗傅里叶分解类似,由平移和调制的高斯窗函数构成小波束分解的框架基本函数.虽然每个小波束具有单一的传播方向,但效率较低.局部余弦基小波束由余弦基调制的钟型窗构成,是正交分解且具有快速算法,并在多个模型上测试了该算法的效率及精度.但局部余弦小波束分解后每个系数都代表以垂直方向为轴的对称传播的波瓣,缺乏单一定义的方向,使其不能用于方向照明等与方向性相关的分析中.为了克服局部余弦基小波束的这一缺陷,毛剑等[7]提出了局部谐和基小波束分解方法.局部谐和基由局部余弦基和局部正弦基线性组合而成,正交且具有单一的传播方向,在复杂模型上比局部余弦基小波束有更高的精度.高静怀等[8]推导了基于一般标架分解的单程波传播算子,并详细讨论了Gabor-Daubechies标架传输算子的渐进展开和高频近似等问题.Ma和Margrave[9]用变窗长Gabor基函数分解波场,使每个窗口内的速度扰动量可以控制,从而提高传播的精度和效率.Pedersen等[10]也提出用横向自适应窗长划分波场,而且根据窗口内速度扰动和精度要求,不同窗口用不同的算子传播.
然而,上述所有局域化方法都是先把波场变换到频率域,沿空间域做局域化分解,然后在频率域偏移成像,只是利用了空间上的局域化.本文探索一种在时间和空间上都局域化的单程波叠前深度偏移方法.物理直观上,时空完全局域化方法会有更多优势.首先,地震波场在时间和空间上都是非平稳信号,时空局域化分解方法更有利于对波场的稀疏表示,也即会对波场有更好的压缩;其次,在沿深度方向的波场延拓过程中,传播一个小的深度后,地震波不会随机地出现在时空平面的任何位置,而完全由波动方程描述的时空关系决定,因此,时空局域化的波场延拓算子应该会更稀疏.近些年引起很多关注的curvelet[11]变换就是一种完全局域化的分解方法,可以用来稀疏表示地震数据和地质结构[12],在理论上被证明有可能是对波的传播算子的最优分解[13].Douma和de Hoop[14]实现了在高频近似下curvelet 域的映射偏移(map migration).Wu等[15-16]利用时间域局部指数标架(称为drumbeat原子)和空间域局部余弦基(beamlet)的张量积,构成了时空完全局域化的dreamlet原子(dreamlet是由drumbeat和beamlet两个词组合而成),并在一些模型上测试了基于该分解的单程波叠前深度偏移.由于局部指数标架[17]是冗余度为2 的紧标架,所以由此构成的dreamlet原子冗余度也为2,在构造波场局域化传播算子时效率难以提高.
本文利用多维局部余弦变换实现dreamlet算法,也即时间和空间上均用正交的局部余弦变换实现局域化分解,构成正交的dreamlet分解原子.正交dreamlet分解不仅可以快速地在时空域和局部时空域转换波场,而且变换域的波场和单程波传播算子均为实数,也即波的传播和成像均在实数域完成,这无疑会加速波的传播和成像过程.本文介绍了正交dreamlet变换及其对时空域波场和频率波数域单程波算子的分解,在二维SEG/EAGE 盐丘模型和SIGSBEE 模型上验证该理论推导的正确性以及在复杂模型和地质结构上的精度.最后是对该方法的一些讨论和结论.
2 正交dreamlet变换及其对时空域波场的分解 2.1 局部余弦基和正交dreamlet原子局部余弦变换是局域化了的离散余弦变换.离散余弦变换[18]是一类线性变换,因为对平滑和相关性较高数据有很强的压缩能力,被广泛的应用于图像编码领域.离散余弦变换具有很多与离散傅里叶变换相似的性质,而且可以利用快速傅里叶变换实现快速运算.离散余弦变换对实数信号的变换输出结果仍为实数信号,对复数信号的实部和虚部分别变换.局部余弦基函数是由平滑、紧支撑的钟函数与一个余弦函数的乘积构成[19],这些基保持了正交性并且具有较小的Heisenberg乘积.
空间域,局部余弦基定义如下:
(1) |
其中,${{\bar x}_n}$为钟函数空间起始位置,${L_n} = {{\bar x}_{n + 1}} - {{\bar x}_n}$为空间长度,m为局部波数指数,Bn(x)是在${{\bar x}_n}$处的钟函数.局部余弦基在波数域的表达式如下:
(2) |
其中,Bn(ξ)为钟函数Bn(x)的傅里叶变换,且珋ξm=π$\left( {m + \frac{1}{2}} \right)/{L_n}$.
在时间域类似,局部余弦基可由时间起始位置珋${{\bar t}_j}$,时间长度Tj=${{\bar t}_{j + 1}}$-${{\bar t}_j}$,和局部频率指数i表示,其时间和频率域表达式分别为
(3) |
(4) |
其中,Bj(ω)是Bj(t)的傅里叶变换,且有$\bar \omega = \pi \left( {{\text{i}} + \frac{1}{2}} \right)/{T_j}$.
时空域dreamlet原子由时间和空间域局部余弦基的张量积构成:
(5) |
dreamlet原子的支集为时间和空间钟型函数的笛卡尔乘积(Cartesianproduct).图 1给出了一种时空域钟型函数的形状.
图 2给出一些正交dreamlet原子在时空域的样本.横轴为空间坐标,dreamlet局域波数由低至高,纵坐标代表时间,局部频率由低至高.这些原子有一定持续时间而且在空间上有一定宽度,并沿空间和时间以一定波数和频率震荡.
反射地震勘探中,地表激发的震源产生能量向地下传播,当遇到地下介质速度或密度变化时,地震波会分裂成两部分,一部分反射到地表被检波器接收,另一部分折射,向更深处传播.地表接收到的地下每一个地层或者散射点的反射信号在地震数据中都分布在一条光滑的曲线上(三维时在时空曲面上).因此,在时空上局域化的分解方法应该会对波场有好的压缩,在传播上又有空间局域化的传播精度.
深度z的时间和空间域波场uz(x,t;z)可以表示为时空局域化dreamlet原子的叠加:
(6) |
其中〈,〉代表内积运算.为了简化表示用${\bar \mu }$=$\left( {{{\bar t}_j},{{\bar \omega }_i},{{\bar x}_n},{{\bar \xi }_m}} \right)$代表局域化指数,${U_{\bar \mu }}\left( z \right)$表示波场的dreamlet分解系数.
我们以二维SEG/EAGE 模型中第175炮数据为例,展示正交dreamlet变换对波场的分解.图 3a为地表检波器记录的地震图.图 3b是正交dreamlet分解系数,时间窗长度为16 个时间采样点(采样间隔为0.008s),空间窗长度为16 个空间采样点(采样间隔为25m).为了显示dreamlet系数分布细节,选择了三个时空窗放大显示.在每个时空窗内,波场的系数主要集中在窗的左上角(低频率低波数),这是由于余弦变换有能量压缩功能.变换后波场系数分布表明,时空局域化的dreamlet分解能够对地震波场实现稀疏表示.在偏移时,只需要传播波场大的dreamlet系数,就可以得到高质量的像.可以有多种方法压缩变换后的波场,但这并非本文的重点,在后面的偏移试验中,我们都以分解后波场系数最大值乘以10-4作为压缩波场的阈值.
不考虑密度变化时,二维标量波波动方程可以表示为
(7) |
其中,u(x,z,t)为压力波场,v(x,z)为波在介质中传播的速度.如果规定深度方向z为波场的延拓方向,而且形式上把对z的二阶偏导看做一阶偏导的平方,式(7)可用平方根公式分解为
(8) |
则全波波动方程可以转化为单程波动方程组:
源波场由方程9a沿深度方向延拓,而地表接收到的波场由方程9b反传至地下,u- 和u+ 分别表示上下行波场.公式中的平方根算子可作为拟微分算子处理,可以由一系列偏微分算子近似展开表示.公式(9)的解为
(10) |
式(10)只是形式地给出了单程波方程的解,该式并没有精确的求解方法,本文通过局部扰动理论求解该方程.局部扰动理论分两步,局部背景速度传播和局部相位屏校正.本文主要推导局部背景下的dreamlet传播算子.
把波场用dreamlet 分解后,可以通过延拓dreamlet 原子来得到更深一层的波场,也即dreamlet原子在介质中的传播必须满足波动方程(10):
(11) |
其中,a(x,t;Δz)表示dreamlet原子dμ(t,x)以其所在处局部匀速背景v(${{\bar x}_n}$,z)下传播深度Δz后的波场.Dreamlet原子在介质中传播一定距离后一般不会再是一个单一的原子,对通过介质演化后的原子波场再分解会得到在新的深度下的dreamlet系数.Dreamlet传播算子是把上一层深度的dreamlet系数映射到当前深度dreamlet系数的权重:
(12) |
式中,${d_{\bar \mu '}}$t,x)表示更深一层的dreamlet系数,传播算子$P_{\bar \mu }^{\bar \mu '}$表示深度z的dreamlet原子${d_{\bar \mu }}$(t,x)至深度z+Δz的dreamlet原子${d_{\bar \mu '}}$(t,x)之间的耦合系数.
利用Parseval定理和傅里叶变换与偏导的对应关系:
(13) |
其中ξ 为水平波数,则dreamlet传播算子
(14) |
real()表示取实部,${d_{\bar \mu }}$(ω,ξ)= bij(ω)bmn(ξ)是dreamlet的频率波数域表达式.从公式推导可以看出,虽然dreamlet传播算子是由频率波数域波的复数传播算子推导得出,但是最后只保留了计算结果的实部.这是因为由两维局部余弦变换构成的正交dreamlet分解的波场系数为实数,在dreamlet域中传播时,待传波场和传播后波场均为实数,那么传播算子也一定是实数.假使传播算子虚部存在值,在传播时也将不起作用,后面的数字试验验证了该推测.把${d_{\bar \mu }}$(ω,ξ)=bij(ω)bmn(ξ)代入(14)得
(15) |
其中$P_{m,n}^{m',n'}\left( {\vartriangle z;\frac{\omega }{{v\left( {{{\bar x}_n},z} \right)}}} \right)$是频率域的小波束(beamlet)传播算子.Dreamlet可以理解为频率域beamlet传播算子增加了时间局域化.
在偏移时,为了提高计算效率,根据整个模型的速度分布,可以提前计算不同参考速度下的dreamlet传播算子,并将其存储于内存中.在每个深度,根据空间窗口的速度在表中挑选与其最接近的算子来传播波场.图 4a显示了在速度为2000 m/s下的dreamlet传播算子,并放大三个时空窗来显示算子的细节分布.从图中可以看出dreamlet传播算子是非常稀疏的,所有的值都集中在两个条带上,而且这些位置遵守波在空间中传播的时间-空间关系.在图 4b和图 4c中显示的是相同速度下从零频至奈奎斯特频率的beamlet传播算子,其中图 4b和4c分别为传播矩阵的实部和虚部.图中明显看出时空都局域化的dreamlet传播算子比仅在空间局域化的beamlet算子更稀疏.如果以传播矩阵中模值(或绝对值)最大值的10-3作为阈值,则dreamlet的稀疏程度(所需存储的系数占整个传播矩阵比例)为7.73%,而相同情况下beamlet传播矩阵稀疏度为35.11%.
炮域叠前偏移先对单炮分别成像,最终的偏移剖面是把各炮偏移结果相叠加.在单炮成像时,炮点处的虚拟源uS(x,z=0,t)与在该点处激发震源地表接收到的地震记录uSR(x,z=0,t)分别延拓(正传与反传)至整个空间.空间中每点的像值由该点处模拟的震源信号和地表上接收到的地震记录反传到该点的信号互相关得到:
(16) |
其中,T是检波器记录的时间长度.相关运算是在比较同一点处地震记录的产生信号和地下结构信息携带信号的相似性,如果在相同时刻(不一定在零时刻)两条记录中都出现子波,则表明该点是一个散射点,而且强度与两者相似程度正相关.
3.4 dreamlet叠前深度偏移流程现将给定模型的dreamlet叠前深度偏移算法(图 5为流程图)总结如下:
(1) 根据模型最大和最小速度选择参考速度间隔,计算dreamlet传播算子并以速度为索引存储于表中,以备所有炮共用;
(2) 对所有炮和所有深度:
(i) 由单炮采集系统分布在dreamlet算子表格中挑选单炮偏移的局部参考速度;
(ii) 由参考速度和真实速度差异实行局部相位屏校正;
(iii) 把时间空间域波场分解至dreamlet域并延拓;
(iv) 变换波场至时间空间域,应用成像条件;
(3) 单炮成像结果叠加至总剖面中.
4 数值算例 4.1 单个dreamlet原子在匀速空间的传播Dreamlet偏移方法中,先对波场进行dreamlet分解,然后用dreamlet传播算子传播dreamlet原子来达到波场延拓的目的,所以单个dreamlet原子在空间中的响应就非常重要,就如同研究空间中点源或者波数域平面波响应一样.图 6a和6c给出了垂直和倾斜向下传播的dreamlet原子在匀速空间三个时刻的响应,作为对比图 6b和6d给出了频率域相同角度beamlet原子的响应.倾斜传播的原子以垂向为对称轴向两个方向传播,这是由于空间局域化的局部余弦基总是有对称的正负波数造成的.从这些响应可以看出,可以把dreamlet看作是许多小的具有一定频带的平面波,在空间和传播方向上都有很好的局域化性质.
近些年,在横向速度有强变化地区的精确成像受到越来越多的关注,如美国的墨西哥湾地区.在该地区,盐体的速度比周围介质高数倍,当地震波穿过这些介质时,盐体像透镜一样把信号能量散开或者汇聚,如何处理这些横向速度变化对盐体下部成像非常关键.二维SEG/EAGE 模型的建立就是为了模拟墨西哥湾地区的地质环境,其速度模型如图 7a.该模型在水平方向上有1200个采样点,深度方向有150个采样点,采样间隔均为25 m.其中最小速度为1524m/s,最大速度为4480 m/s,盐丘速度大约是周围介质速度3倍.偏移时,该模型的传播算子用了25个参考速度,从最小速度至最大速度等间隔分布.图 7b为dreamlet叠前深度偏移结果.我们看出该模型的主要结构,如岩体不规则的上边界,模型底部的水平基线,以及盐下的断层结构成像都准确、清晰而且合理.
在用单程波延拓波场偏移成像过程中,用于成像浅层结构的地震信号对深层结构不再有贡献.Dreamlet波场延拓的速度和有效信号量成正比,因此为了提高偏移速度,我们以速度模型中的最大和最小速度估算走时,去掉检波器波场中用来成像浅层结构的信号,而不影响成像质量.图 8 输出了第175炮检波器波场在不同深度的地震图.从图中看出,随着深度的增加,地震记录中的有效信号时间长度在不断减小,这就意味着随着深度增加,dreamlet波场延拓效率会提高,这是时间局域化的一个显著优势.
SIGSBEE2A 模型横向速度变化剧烈,盐丘下方有许多薄的反射层,而且,速度模型中有两排小的散射点用以测试偏移算法的聚焦能力,用来偏移成像的速度模型如图 9a所示.该模型水平方向有2133个采样点,深度方向有1200个采样点,采样间隔分别为11.43m 和7.62 m.每炮有384个检波器,检波器时间采样为0.008s,时间长度为12s.图 9b 是dreamlet以速度9a偏移成像的结果,其中平反射层,断层以及盐丘的边界和底部的水平层的结构都成像清楚.另外,两排散射点也成像清晰且在正确的位置.该模型偏移时,dreamlet背景传播算子有50个参考速度,同样从最小至最大速度均匀分布.作为对比,图 9c给出了beamlet偏移结果.Dreamlet偏移是在频率域beamlet偏移的基础上增加了时间局域化,而二者在空间上的分解完全相同,因此两种偏移方法具有相同的适应速度横向变化的能力.在SIGSBEE2A合成数据模型测试算例中,两种算法的得到了类似的偏移结果.
地表接收到的关于地下不连续层或者散射点的信号在地震数据中分布在一个时空光滑平面上(二维时为曲线).由于子波持续期相对于整个地震记录时间很短,对时间轴实现时间频率局域化分解能够更有效地表示反射信号.偏移成像中,最理想的变换或分解方法要同时能对地震数据、像结构以及偏移算子均有好的稀疏表示效果.基于局部扰动理论,本文利用局部余弦基的张量积构成在时间和空间上完全局域化的正交dreamlet原子,能够分别对地震数据和单程波传播算子进行稀疏表示.而且分解后的波场和传播算子都为实数,也即偏移和成像过程均在实数域完成.然而本算法中关于局部速度扰动校正部分仍没有很好的解决,在每个深度,都把波场变换到频率域后采用与文献[6]中相同的方法进行局部相位屏校正,然后再将校正后的波场变换回时间域,因此占用了大量的计算时间.在偏移过程中,虽然检波器波场的有效时间长度随着深度的增加在不断减小,但是源波场在传播过程中随着深度的增加会在空间上不断扩散,从而增加计算量.基于这些特点,后续研究中拟将dreamlet算子与观测系统沉降观念相结合,进一步提升偏移效率.另外,对变换域波场和传播算子采用何种阈值算法压缩也是一个重要的课题,这些都将是下一步研究的重点.
致谢非常感谢谢小碧、贾晓峰、曹军、闫蕊、耿瑜、毛剑、陈文超关于本项工作的讨论和建议,感谢评审专家对本文提出的宝贵修改意见.
[1] | Steinberg B Z. Evolution of local spectra in smoothly varying nonhomogeneous environments local canonization and marching algorithms. J. Acoust. Soc. Am. , 1993, 93(5): 2566-2580. |
[2] | Steinberg B Z, Birman R. Phase-space marching algorithm in the presence of a planar wave velocity discontinuity-A qualitative study. J. Acoust. Soc. Am , 1995, 98(1): 484-494. |
[3] | Wu R S, Wang Y Z, Gao J H. Beamlet migration based on local perturbation theory. // 70th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 2000: 1008-1011. |
[4] | Wu R S, Chen L. Wave propagation and imaging using Gabor-Daubechies beamlets. // Theory and Computational Acoustics, World Scientific, New Jersey, 2002: 661-670. |
[5] | Chen L, Wu R S, Chen Y. Target-oriented beamlet migration based on Gabor-Daubechies frame decomposition. Geophysics , 2006, 71(2): s37-s52.. DOI:10.1190/1.2187781 |
[6] | Wu R S, Wang Y Z, Luo M Q. Beamlet migration using local cosine basis. Geophysics , 2008, 73(5): 207-217. DOI:10.1190/1.2969776 |
[7] | 毛剑, 吴如山, 高静怀. 利用局部谐和基小波束的高精度叠前深度域偏移成像方法研究. 地球物理学报 , 2010, 53(10): 2442–2451. Mao J, Wu R S, Gao J H. High-accuracy prestack depth migration and imaging using local harmonic basis beamlet. Chinese J. Geophys. (in Chinese) , 2010, 53(10): 2442-2451. |
[8] | 高静怀, 周艳辉, 毛剑, 等. 局部角度域波传播步进算法研究. 地球物理学报 , 2007, 50(1): 248–259. Gao J H, Zhou Y H, Mao J, et al. A wave propagation method in local angle domain. Chinese J. Geophys. (in Chinese) , 2007, 50(1): 248-259. |
[9] | Ma Y W, Margrave G F. Seismic depth migration with the Gabor transform. Geophysics , 2008, 73(3): s91-s97. DOI:10.1190/1.2903821 |
[10] | Pedersen Ø, Brandsberg-Dahl S, Ursin B. Seismic imaging using lateral adaptive windows. Geophysics , 2010, 75(2): s73-s79. DOI:10.1190/1.3359799 |
[11] | Candès E J, Demanet L, Donoho D, Li Y. Fast discrete curvelet transforms. SIAM Multiscale Modeling and Simulation , 2006, 5(3): 861-899. DOI:10.1137/05064182X |
[12] | Herrmann F. Optimal seismic imaging with curvelets. // 73th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 2003: 997-1000. |
[13] | Candès E J, Demanet L. The curvelet representation of wave propagators is optimally sparse. Communications on Pure and Applied Mathematics , 2004, 58: 1472-1528. |
[14] | Douma H, de Hoop M V. Leading-order seismic imaging using curvelets. Geophysics , 2007, 72(6): s231-s248.. DOI:10.1190/1.2785047 |
[15] | Wu R S, Wu B Y, Geng Y. Seismic wave propagation and imaging using time-space wavelets. // 78th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 2008: 2983-2987. |
[16] | Wu R S, Wu B Y, Geng Y. Imaging in compressed domain using dreamlets. // CPS/SEG Beijing' 2009, International Geophysical Conference, Expanded Abstracts, ID: 57. |
[17] | Mao J, Wu R S, Gao J H. Directional illumination analysis using the local exponential frame. Geophysics , 2010, 75(4): s163-s174. DOI:10.1190/1.3454361 |
[18] | Püschel M, Moura J. The algebraic approach to the discrete cosine and sine transform and their fast algorithms. SIAM Journal of Computing , 2003, 32(5): 1280-1316. DOI:10.1137/S009753970139272X |
[19] | Coifman R R, Meyer Y. Remaxques sur 1'analyse de Fourier a fengtre. Comptes Rendus de l' Acad6mie des Sciences , 1991, 312: 259-261. |