2016, Vol. 59
2. 海洋国家实验室海洋矿产资源评价与探测技术功能实验室, 青岛 266071
2. Laboratory for Marine Mineral Resources, Qingdao National Laboratory for Marine Science and Technology, Qingdao 266071, China
地下强反射界面产生的多次波一直是困扰地震勘探数据处理以及解释的棘手难题.为了消除多次波对一次反射波成像的影响，常规数据处理中通常将其视为噪声，并根据它与一次波的时差差异以及周期性基于滤波理论发展了一系列多次波去除方法(Taner and Koehler，1969；Ryu，1982；Hampson，1986；等).随着研究的深入，基于波动理论的多次波预测相减法(Wiggins，1988；Berkhout and Verschuur，1997； Verschuur and Berkhout，1997；Weglein et al.，1997；等)为有效去除多次波开辟了一条新的途径，其中，Berkhout和Verschuur(1997)提出的基于反馈环理论的表层多次波去除技术(SRME)在目前的工业生产中有着最为广泛的应用.
然而，考虑到多次波也是来自地下界面的真实反射，且相比一次波具有地下传播路径长、覆盖区域广及反射角度小的特点，越来越多的地球物理学家开始将研究的目光投向多次波的有效利用中，并通过多次波成像技术研究挖掘多次波中蕴含的地下界面有效信息，获取更广的照明范围和更多的地下构造信息.Berkhout和Verschuur(1994)最先提出将表层多次波作为二次震源进行偏移的方法；Guitton(2002)用一次波和多次波记录作为面炮震源代替震源子波，将多次波作为接收记录，采用一次波单程波偏移方法对多次波进行成像；Liu等(2011)结合逆时偏移基于同样的方法对预测的多次波进行成像；Shan(2003)提出多次波广义炮检偏移概念，通过将一次波和对应的表层多次波进行互相关构建准一次波，然后对准一次波采用常规的炮检偏移间接实现多次波的成像，由于构建的准一次波与一次波形态相似，因此可将其用于补充由于常规观测系统限制造成的远或近炮检距数据的缺失(Shan and Guitton，2004)或者成像区域的缺失(郭书娟等，2012)；Berkhout和Verschuur(2003)提出“聚焦变换(focal transform)”这一概念，并在聚焦变换域内实现多次波的成像(Berkhout and Verschuur，2004，2006)，但是其研究仅限于多次波成像，而未对多次波的其他应用进行深入分析.也有很多学者尝试将多次波和一次波同时进行偏移成像，如He和Schuster(2003)，Brown和Guitton(2005)，郭书娟等(2011)，Berkhout和Verschuur(2011)等.在多次波发育的地震资料中，尤其是海上勘探资料，往往表层多次波发育且能量较强，而层间多次波能量较弱且不发育，因而上述方法研究都存在如下假设：将层间多次波视为“一次波”，仅针对表层多次波的有效利用开展研究，本文也是基于这一假设展开.
通过多次波成像方法研究可以获取地下一次波照明无法到达的地下构造信息.但是，多次波的有效利用并不局限于此，不同于上述思路，本文根据多次波与一次波传播路径关系，结合Berkhout和Verschuur(2003)提出的聚焦变换，在聚焦变换域内利用多次波提高地震数据的分辨率，以崭新的角度实现了对多次波有效应用的深入挖掘.该方法对于地震勘探尤其是海上地震勘探中受震源子波及其他采集因素影响造成的低分辨率地震数据的处理解释具有十分重要的研究意义.
本文首先通过公式推导得出聚焦变换域内多次波表现为原始数据的多维子波反褶积这一关键性结论，从理论上证明了本文方法利用多次波可以有效提高数据的分辨率.然后根据聚焦域内一次波和多次波能量的分布特点，在τ-p域引入Fomel(2009)提出的整形正则化非稳态回归自适应匹配滤波方法来进行聚焦域内高分辨率数据的提取.采用一次波处理方法对提取的高分辨率数据进行处理即可获得高分辨率的处理结果.最后通过模型试算以及实际资料处理验证本文方法的有效性、适应性和实用性.
2 聚焦变换提高分辨率原理2003年，Berkhout提出了聚焦变换这一概念，并给出了聚焦变换正反变换的表达式如下：
|
(1a) |
|
(1b) |
其中，矩阵 P 表示数据体矩阵，矩阵 G 表示传播算子矩阵， F = G －1，聚焦变换可定义为每个空间频率分量上矩阵的相乘.根据聚焦域能量分布特点(图 1所示)，一系列相关技术应运而生.聚焦变换的提出为去噪、数据规则化、信号重构以及多次波压制等处理过程的实现提供了一个新的途径.
|
图 1 聚焦变换能量分布图(Berkhout and Verschuur，2004) Fig. 1 Energy distribution of focal transform(Berkhout and Verschuur，2004) |
聚焦变换的提出最先受启发于数据驱动的表层多次波去除(SRME)技术，因此，为了深入推导含多次波数据的聚焦变换表示，获取聚焦域的高分辨率地震数据，本文首先对SRME理论进行简要的回顾.基于反馈模型推导的SRME核心公式可采用频域数据矩阵表示如下，其中数据矩阵的求取可参考Berkhout(1982)：
|
(2) |
其中， P 表示含多次波的地震数据矩阵，Δ P 表示一次波矩阵，M 表示多次波矩阵，A 表示地表算子.由 M =Δ PAP 可知，多次波可由原始数据矩阵 P 乘以算子矩阵Δ PA 得到.
对(2)式进行推导，并结合级数展开可得
|
(3) |
其中，Δ P 表示一次波(即零阶多次波)，(Δ PA)Δ P 表示一阶多次波，(Δ PA)2Δ P 表示二阶多次波，以此类推.由(3)式可知，对一次波数据矩阵Δ P 乘以加权的一次波算子(Δ PA)即可将一次波转化为一阶多次波，同样对一阶多次波数据矩阵乘以该加权算子可实现一阶多次波向二阶多次波的转化，以此类推.将(3)式代入多次波的表达式中，得
|
(4) |
上述过程中与算子(Δ PA)相乘是实现多次波升阶的过程.与此相反，利用一次波和多次波传播路径上的关系，通过引入聚焦变换也可以实现多次波的降阶过程.定义聚焦变换算子如下：
|
(5a) |
|
(5b) |
此时，正聚焦变换表示将由一次波构建的正聚焦变换算子 F 与含多次波的原始地震数据矩阵 P 进行相乘，变换到聚焦域.而反聚焦变换则表示将聚焦域数据变换回原始数据域.
将式(5b)和式(3)代入式(1a)，得
|
(6) |
其中， I =Δ P －1Δ P，表示聚焦域内位于原点周围的聚焦点.由于 A 为与地表有关的加权因子，因此(6)式可以看作是多次波降阶的过程，即一阶多次波转化为一次波，二阶多次波转化为一阶多次波，以此类推.
将式(2)代入式(6)，可得
|
(7) |
式中第一项为聚焦域内一次波产生的聚焦点，而第二项则为聚焦域内由多次波构建的新数据，记为 P ′.下面将结合反褶积理论以及Berkhout的反馈模型对第二项展开深入讨论，给出本文利用多次波提高地震数据分辨率的理论证明.
野外采集中受复杂子波作用以及干扰影响，采集的实际地震记录往往分辨能力较低，在地质界面上各反射波互相叠加、彼此干涉，所以无法获取准确的地质界面.为了解决这一问题，反褶积被广泛地用于提高地震数据的分辨率中(李振春和张军华，2004).
在不考虑干扰的情况下，地震记录x(t)可以表示为地震子波b(t)与反射系数ξ(t)的褶积，
|
(8) |
将其转化到频率域，有
|
(9) |
令A(ω)=1/B(ω)，则可得反褶积的频域表示形式如下：
|
(10) |
其中，A(ω)为反子波的频域表示，在时间域代表子波反褶积滤波器.(10)式即为传统反褶积模型，基于反射系数满足随机白噪这一假设，通过求取反子波实现地震数据的高分辨率转换，但是实际上反射系数并不符合白噪特征，于是很多学者基于上述理论发展了多种改进反褶积模型，尽管如此，各种改进方法的应用范围仍存在一定的局限性.而本文方法利用多次波提高地震数据分辨率则是基于波动理论模型，相比传统的反褶积模型具有更广的适应范围.
根据Berkhout(1982)由波动理论推导的反馈模型，有
|
(11a) |
|
(11b) |
其中，Δ P －表示上行一次波，X 为地下传播算子矩阵，S +为下行的震源波场，它的每一列代表一个震源向量，M －表示上行的多次波，R －为自由地表反射系数矩阵，P －为上行总波场.
而上行多次波又可以表示为
|
(12) |
同式(2)，A 为地表算子，将式(11a)和(11b)代入式(12)，推导可得
|
(13) |
通常，可以将自由地表反射系数矩阵 R －表示为负的单位矩阵- I，式(13)可以重写为
|
(14) |
由上述理论推导可知，地表算子 A 即为负的逆震源波场的频域表示，其每一列表示为负的反子波频域形式，在时间域内将表示子波反褶积滤波器.对比式(10)可知，经聚焦变换后式(7)中第二项由多次波恢复的数据相当于对原始数据进行了子波反褶积，这即为聚焦域内利用多次波可以提高地震数据分辨率的原因所在.与常规反褶积单道处理不同，该方法输出的每一道结果中所有的炮记录都参与了运算，因此该过程相当于一个多维子波反褶积过程，在提高垂向分辨率的同时在空间上加以约束，有助于提高数据的横向分辨率.这对于低分辨率的地震数据处理是极为有利的.需要指出，为了与传统反褶积理论保持一致，在得到高分辨率数据以后还需要通过线性运算进行极性校正，去除负号的影响.
相比传统反褶积模型，本文基于Berkhout的波动方程模型推导得出的多维反褶积理论不存在假设、不需要求取子波，且可以适用于任何复杂情况，其适用性更广.
3 整形正则化非稳态回归自适应匹配滤波提取高分辨率数据由于聚焦域内一次波的能量并不能完全聚焦在一个点上，而是聚焦在时间t=0、偏移距x=0处附近，表现为两条交叉的射线.在地下复杂以及浅水情况下，一次波聚焦能量与多次波构建的高分辨率数据可能会产生能量干涉；同时聚焦变换过程中受反问题求解影响，在聚焦域内会引入大量噪声，因此，简单的切除显然无法实现高分辨率数据的合理提取.如何提取出由多次波构建的高分辨率数据成为重中之重.
考虑到聚焦域内能量聚焦的特点，可选择将聚焦域结果先映射到τ-p域，然后通过自适应匹配滤波方法提取出高分辨率数据.τ-p变换可将与角度有关的一次波信息置于τ=0处(de Bruin et al.，1990)，这对于聚焦域内存在一次波和多次波能量干涉数据的自适应匹配是极为有利的；此外，自适应匹配滤波要求数据具有较高的信噪比，τ-p域转换可以有效压制聚焦域内产生的噪声，提高数据的信噪比，进而提高匹配的精度.为了尽可能减少τ-p正反变换引入的误差，本文采用高分辨率τ-p变换实现数据映射.前人研究表明高分辨率τ-p变换在接受精度范围内可认为是保幅的，所以由此产生的误差可以忽略.综合考虑以上因素可知，τ-p域映射可以很好地辅助后续的自适应匹配，利于高分辨率数据的提取.
目前发展的自适应匹配滤波方法大都基于信号平稳假设，但是地震信号通常无法满足这一假设.因此，在进行匹配滤波时需要选取时窗，而时窗大小控制不当极易产生不稳定、有效能量损失或者噪声残留等问题，造成波场分离精度受限，无法对子波进行有效校正，因而不能保证数据的最佳匹配，这对于高分辨率数据的提取是极为不利的(Abma et al.，2005).
于是，本文采用非稳态回归滤波方法来进行聚焦域内波场的分离，并引入Fomel(2009)提出的整形正则化条件.该方法通过设计关于时间和空间变化的非稳态匹配滤波算子，不需要将数据分时窗处理，匹配滤波算子可以适应地震数据时间和空间的变化；通过整形正则化对滤波系数的连续性和平衡性进行约束，可以避免不稳定问题的产生，同时也提高了计算效率.
对于本文中聚焦域内的波场分离，基于稳态回归的自适应匹配滤波算法的误差函数可以表示为
|
(15) |
其中，x 为多维空间坐标，p′(x)为聚焦域数据，pk(x)是由原始地震数据经不同空间移动后得到的序列，ak为滤波系数.由式(15)可看出滤波系数并不随着信号坐标的变化而改变.为了尽可能满足稳态假设，在处理中需要分时窗处理，尽管分时窗处理对结果有一定改善，但是却不能从根本上解决问题.
而非稳态回归自适应匹配滤波算法则考虑了滤波系数的空间变化，其误差函数可表示为
|
(16) |
由于滤波系数随坐标变化，因而未知变量多于方程个数，(16)式的求解是病态的，为了解决这一问题，必须对滤波系数进行附加约束，这就需要加入正则化.本文引入了Fomel(2009)提出的整形正则化方法，其对应的目标函数表示如下：
|
(17) |
其中，S 为正则化算子，λ为正则化系数，该方法中可将正则化问题看作是模型域的平滑操作.相比常用的Tikhonov正则化(Tikhonov，1963； Engl et al.，1996)，该方法更利于正则化算子和正则化系数的选取，且反演速度快，为解决不适定的反问题提供了更为简便快速的方法.
采用整形正则化非稳态自适应匹配滤波后，即可将聚焦点(t=0)周围的能量与其他的反射能量区分开，提取出由多次波构建的高分辨率地震数据，对该数据进行多次波去除以后再进行叠加或者成像处理，即可获取高分辨率的处理结果.
|
图 2 高分辨率数据提取示意图 Fig. 2 Illustration of extraction of data with high resolution |
图 2给出了本文方法实现的流程图.首先通过SRME提取一次波，然后由一次波构建正聚焦变换算子，并对原始数据的频域数据矩阵进行正聚焦变换获取聚焦域内结果 Q，最后在τ-p域通过本文引入的整形正则化非稳态自适应匹配滤波即可获取高分辨率数据 P ′，反τ-p变换得到利用多次波获取的高分辨率数据.需要指出，为了获取高精度的处理结果，聚焦变换通常需要进行多次迭代.采用SRME可以获取较精确的一次波，经1~2次聚焦变换的迭代过程即可获取较高精度的聚焦变换结果.随着计算机技术的不断发展，SRME的实现过程已不再受计算量的限制，但是在SRME较难实现的情况下，也可以通过其他方法来获取一次波，如滤波法，只是相比SRME，一次波估计精度较低，此时需要通过多次迭代来获取高精度的聚焦变换结果.
4 数值试算为了对本文方法的正确性进行合理的验证，文中分别对简单平层模型、复杂模型以及实际资料进行了数值试算.通过简单平层模型试算对本文方法的实现过程及优势体现展开详细的分析讨论，验证方法的有效性；通过复杂模型(含复杂构造、特殊储集体及薄层构造)试算对该方法的适应性进行深入讨论；最后通过实际资料试处理验证本文方法的实用性.
4.1 平层模型试算首先，采用简单的平层模型来验证本文方法的有效性.速度模型如图 3a所示，该模型包含2个相邻的反射界面，横纵向均为101个采样点，空间采样率为10 m，从上到下三层介质速度分别为1500 m·s-1，1750 m·s-1，2100 m·s-1.假设地表为自由反射界面，采用有限差分方法进行正演模拟，在地表以10 m 间隔均匀分布101个炮点，以10 m间隔均匀分布101道接收.为了有效验证本文方法在提高分辨率方面的优势，选取较复杂的子波如图 3b所示，该子波为2014年SEG网站公开数据库提供的用于2D海上勘探数据模拟所采用的远场子波，子波采样率为0.666 ms，由图可见子波的旁瓣较为复杂，这会严重影响地震数据的分辨率，进而影响处理结果的质量.
|
图 3 平层模型速度场(a)及正演子波(b) Fig. 3 Flat layer model(a)velocity model and(b)wavelet |
图 4展示了正演炮记录聚焦变换的结果.图 4a为正演的单炮记录，时间采样率为8 ms.从图中可以看出，由于子波复杂，炮记录中同相轴的分辨率较低，两界面对应的一次反射波能量互相干涉，这对于后续的成像及解释工作是极为不利的.同时由于自由地表的存在，在1.1~1.4 s之间产生了明显的多次波反射.首先对该记录采用SRME提取一次反射波，结果如图 4b所示，将其作为聚焦算子，根据式(1a)进行聚焦变换并经极性校正后得到聚焦域结果如图 4c所示，由图可见，在t=0及其附近聚焦有较强的一次反射能量，而多次波则在下部呈现为类似一次波形态的新记录，且可以看出新构建的记录中同相轴相比图 4a中明显变细，但是一次波聚焦能量并不是聚焦在一点而是呈放射性分布，且在聚焦域内引入了部分噪声，这些噪声考虑是在聚焦变换过程中由于矩阵求逆引起的.如何去除一次波能量和噪声，有效提取新构建的高分辨率数据给我们提出了新的问题.根据一次波和多次波在聚焦域的表现形式，考虑先将聚焦域结果通过高分辨率τ-p变换映射到τ-p变换域，如图 4d所示，τ-p域内一次波分布在τ=0处，多次波则表现为不同曲率的曲线.
|
图 4 平层模型聚焦变换结果 (a)单炮记录；(b)聚焦算子；(c)聚焦域结果；(d)图(c)数据τ-p域结果. Fig. 4 The focal transfrom result of flat layer model (a)Single shot record;(b)Primaries as the focal operator;(c)Result in the focal domain;(d)Result in the τ-p domain of figure(c). |
为了在有效保存数据分辨率的前提下提取出聚焦变换域内多次波构建的新数据，本文采用非稳态回归自适应匹配滤波方法，并结合Fomel(2009)的整形正则化算法在τ-p域内进行聚焦域数据与原始数据的自适应相减，提取出利用多次波构建的高分辨率数据.为了验证该方法的优势，本文分别采用常规基于稳态回归的自适应匹配滤波方法和本文方法对模型数据进行了试算，结果如图 5所示.图 5a为原始单炮记录高分辨率τ-p变换结果，受子波旁瓣影响，两界面对应的一次波反射能量干涉在一起，分辨率较低，两同相轴较难分辨.图 5b则为该炮记录对应聚焦变换结果的τ-p域映射，由图可见，聚焦变换产生的噪声在τ-p变换后得到了有效的压制，这对于获取高精度的匹配滤波结果是十分有利的.将图 5b分别采用常规的稳态回归自适应匹配滤波和本文引入的整形正则化非稳态回归自适应匹配滤波方法与原始数据进行匹配相减，提取出利用多次波构建的高分辨率数据，结果如图 5c和图 5d所示.对比两图可见，常规的稳态回归自适应匹配滤波和本文方法都可以将多次波构建的新数据与一次波聚焦能量分离开来，但是常规方法提取的多次波不能有效保存聚焦域数据的高分辨率，而本文方法则可以对子波进行较好的校正，有效保存数据的分辨率.图 5e和图 5f分别为采用上述两种方法提取高分辨率数据后的剩余能量，可明显看出本文方法相比常规方法可以更好地实现高分辨率数据与一次波聚焦能量的分离.将图 5c和5d结果进行反τ-p变换返回原数据域，所得结果如图 5g和5h.对比两图可明显看出，采用本文方法提取的由多次波构建的新数据相比常规匹配滤波方法提取的数据具有更高的分辨率，一次波反射轴得到了更好的刻画，同相轴变细，不再存在能量的相互干涉，且在聚焦变换域内产生的噪声得到了有效的去除.上述结果对本文方法的有效性进行了验证，并证明了如下结论：基于聚焦变换利用多次波可以构建出高分辨率的地震数据，而采用整形正则化非稳态回归自适应匹配滤波方法可以有效提取出聚焦域内由多次波构建的高分辨率数据.
|
图 5 常规自适应匹配滤波与本文整形正则化非稳态回归自适应匹配滤波结果对比 (a)单炮记录τ-p域结果；(b)聚焦变换结果τ-p域映射；(c)常规自适应匹配滤波结果；(d)正则化非稳态回归匹配结果； (e)常规匹配后误差；(f)本文方法匹配后误差；(g)常规匹配后多次波构建的记录；(f)本文方法匹配后多次波构建的记录. Fig. 5 Result comparison between conventional adptive match filtering and non-stationary regression adptive match filtering with sharping regularization (a)Shot record after τ-p transform;(b)Focal domain data after τ-p transform;(c)Result after conventional adptive match filtering;(d)Result after non-stationary regression adptive match filtering with sharping regularization;(e)Error after filtering by conventional filtering;(f)Error after filtering by method in this paper;(g)Record transformed from the multiples by conventional filtering;(f)Record transformed from the multiples by method in this paper. |
为了更好地突出本文方法的优势，下面将进一步对处理结果进行对比分析.对原始数据和新提取的高分辨数据分别进行叠加和偏移，结果如图 6所示.图 6a和6b分别为对原始数据和本文利用多次波提取的高分辨率数据进行多次波去除、动校正、叠加处理后得到的叠加剖面.由图可见，原始数据的叠加剖面上两界面对应的同相轴受子波影响，分辨率较低，而采用本文方法利用多次波提取的高分辨率数据对应的叠加剖面上同相轴明显变细，分辨率得到了有效的提高.对比两叠加剖面的频谱图可见，采用本文方法获取数据的叠加结果中低频成分和高频成分都得到了较好的补充，频谱明显变宽，分辨率提高得到了有效的印证.同样，对比图 6c和6d的逆时偏移(reverse time migration，RTM)结果，可得到同样的结论，图 6c中受子波影响，偏移结果中两界面较难分辨，这可能会造成错误的解释结果，而图 6d中两界面则较为清晰，分辨率较高.
|
图 6 原始数据与本文聚焦域提取数据叠加和偏移结果对比 (a)原始数据叠加结果；(b)聚焦域提取数据叠加结果；(c)原始数据RTM结果；(d)聚焦域提取数据RTM结果. Fig. 6 Comparison of stack result and migration result between original data and separated data form the focal domain by the method in this paper (a)Stack section of the original data;(b)Stack section of the data separated by the method in this paper; (c)RTM result of the original data;(b)RTM result of the data separated by the method in this paper. |
采用本文方法对含特殊构造的复杂模型进行高分辨率数据提取，并重点对复杂构造、特殊储集体和薄层的成像结果进行分析讨论，进一步验证本文方法的适应性.
选取模型速度场如图 7a所示.该模型中第一个界面用于模拟海底，模型中包含背斜构造、上覆层状构造、断层构造和深部的薄层，还有裂缝和孔洞构造，横向301个采样点，纵向251个采样点，空间采样率为10 m，介质速度范围为2400~5500 m·s-1.假设地表为自由反射界面，采用图 3b所示子波进行正演模拟，图 7b为正演的炮记录，而图 7c则为采用本文方法利用多次波在聚焦域提取的新的炮记录，对比图 7b和图 7c，可明显看出利用多次波构建的炮记录中同相轴更细，具有更高的分辨率.由于波场较为复杂，为了更好地说明这一问题，对正演数据进行偏移，从偏移结果上来进行更好的对比分析.首先对原始数据以及本文方法提取的数据进行多次波压制，然后再采用RTM进行偏移，结果分别如图 8a和8b所示，对比两图可看出图 8b中界面同相轴分辨率明显高于图 8a，海底界面、背斜顶部和侧翼以及上覆构造的同相轴分辨率都得到了明显的提高，尤其是在深度2 km左右薄层构造对应的两个界面同相轴，由于原始数据受子波影响，成像结果中两同相轴能量干涉在一起，无法进行有效区分，而图 8b采用本文方法提取的高分辨率数据对应的偏移结果则能够较好地区分薄层构造.此外，从背斜侧翼的细节构造中可以看出成像结果的横向分辨率也得到了较好的改善.图 8c和图 8d分别对该模型中特殊储集体成像结果进行进一步对比，由图可见，图 8d中方框内的孔洞构造、裂缝构造以及断层构造都较图 8c得到了更好的刻画，同相轴更细，能量更为收敛，分辨率得到了一定的提高.需要指出，考虑模型的复杂程度，在进行自适应匹配过程中，需要选取合适的pk(x)序列长度，否则容易造成提取数据的分辨率受损或者在成像结果中引入噪声.
|
图 7 复杂模型(a)速度场，(b)原始炮记录，(c)聚焦变换域多次波转化炮记录 Fig. 7 Complex model(a)velocity model，(b)original shot record，(c)shot record transformed from multiples in the focal domain |
|
图 8 复杂模型偏移结果 (a)原始炮记录偏移结果；(b)本文方法聚焦域提取记录偏移结果；(c)图(a)局部放大图；(d)图(b)局部放大图. Fig. 8 Migration result of complex model (a)Migration result of original shot record;(b)Migration result of transformed shot record in focal domain by the method in this paper; (c)Enlarged partial view of figure(a);(d)Enlarged partial view of figure(b). |
最后将本文方法应用于实际资料，进一步验证该方法的实用性.对多次波发育的某海上实际资料采用本文方法提取高分辨率数据并将其与原始数据进行对比.文中为了更直观地对比原始数据跟本文方法提取数据的分辨率，分别抽取了原始炮记录和利用多次波构建记录的近偏移距剖面，并对其进行对比分析，结果如图 9所示.图 9a为原始数据x= 200 m的共偏移距剖面，从图中可看出在3.2~3.8 s 之间存在发育的多次波，图 9b为采用本文方法提取高分辨率数据并进行多次波压制以后提取的x=200 m 的共偏移距剖面，从图中可看出多次波得到了有效的压制，且整体剖面的信噪比得到了一定的提升，这是由于聚焦变换中信噪分布区域不同以及τ-p域映射作用，在提取高分辨率数据的同时，原始数据中的噪声也得到了一定的压制.对图 9c和9d中局部放大图的细节构造进行对比，可明显看出上部方框内的海底界面同相轴以及在时间t=2.4 s处方框内界面同相轴的复杂旁瓣都得到了有效的压缩，分辨率得到了明显提高.而对比2~2.2 s处方框内的同相轴可明显看出，图 9d中不仅垂向分辨率得到了一定的提高，同时，同相轴在空间上的连续性也得到了一定的加强，这是由于本文方法在进行聚焦变换的过程中每一炮都参与了运算，相当于对空间进行约束，在一定程度上可以提高数据的横向分辨率.
|
图 9 实际资料近偏移距剖面对比图 (a)原始数据x=200 m共偏移距剖面；(b)本文方法提取数据x=200 m共偏移距剖面；(c)图(a)局部放大图；(d)图(b)局部放大图. Fig. 9 Comparison of near offset section of real field data (a)x=200 m common offset section of original data;(b)x=200 m common offset section of data transformed by the method in this paper; (c)Enlarged partial view of figure(a);(d)Enlarged partial view of figure(b). |
目前对多次波的有效利用都是围绕多次波成像技术研究开展的，利用多次波成像来获取更广的地下照明和更多的地下信息.不同于这一常规思路，本文从利用多次波提高分辨率这一角度对多次波有效利用进行了深入挖掘.研究发现，聚焦变换可以实现多次波的降阶过程，即在聚焦域内可以由多次波恢复出原始数据的形态，通过理论推导得出多次波在聚焦域内表现为原始数据的多维子波反褶积这一重要结论，从理论上对该思路的可行性予以证明，同时还发现，该多维子波反褶积过程中每一炮都参与了运算，增加了空间约束，因此，采用本文方法利用多次波不但可以提高数据的垂向分辨率，而且对于横向分辨率的提高也有一定的辅助作用.相比常规的反褶积模型，该方法基于波动理论推导实现，不存在假设且不需要求取子波，可以适用于任意复杂情况.
为了有效提取由多次波构建的高分辨率数据，解决聚焦域内一次波聚焦能量和噪声的去除问题，文中首先将聚焦域结果采用高分辨率τ-p变换映射到τ-p域，然后采用引入整形正则化的非稳态回归自适应匹配滤波方法进行高分辨率数据的提取.不同于常规基于稳态回归的自适应匹配滤波方法，文中方法匹配因子随空间变化，避免了由于时窗选取带来的一系列问题.通过整形正则化进行反演问题约束，为解决不适定的反问题提供了更为简便快速的方法.采用该方法可以较好地实现高分辨率数据的提取，有效保存数据的分辨率，同时可以压制聚焦域由聚焦变换产生的噪声.
最后通过简单平层模型，复杂模型以及实际资料试算对本文方法的有效性、适应性和实用性进行了验证.处理结果显示，对于多次波发育的地震资料，本文方法可以有效消除复杂子波的影响，结合本文引入的非稳态回归自适应匹配滤波可以在聚焦域内提取出高分辨率地震数据，这对于后续地震资料的处理和解释有着十分重要的意义.
利用多次波不仅可以获取更多的地下照明信息和构造信息，还可以通过本文方法利用多次波提高地震数据的分辨率，如何将这两种优势结合是值得思考的一个问题.此外，本文方法实现涉及到一次波提取、聚焦变换、高分辨率τ-p变换、非稳态回归自适应匹配滤波等计算过程，尽管随着计算机技术的不断发展，这些技术的实现变得切实可行，但是仍需在算法优化以及合理运用并行技术等方面开展进一步研究，尽可能的节省计算成本.
致谢感谢开源软件Madagascar提供的绘图支持和软件支持，感谢SEG网站提供的数据支持，特别感谢两位评审专家为本文完善提出的宝贵意见，感谢贵刊编辑及王鹏博士对本文修改给予的帮助
| Abma R, Kabir N, Matson H K, et al. 2005. Comparisons of adaptive subtraction methods for multiple attenuation. The Leading Edge , 24(3): 277–280. | |
| Berkhout A J. Seismic Migration: Imaging of Acoustic Energy by Wave Field Extrapolation, A: Theoretical Aspects. New York: Elsevier, 1982 : 211 -218. | |
| Berkhout A J, Verschuur D J. 1994. Multiple technology: Part 2, migration of multiple reflections.//64th Annual International Meeting, SEG, Expanded Abstracts, 1497-1500. | |
| Berkhout A J, Verschuur D J. 1997. Estimation of multiple scattering by iterative inversion, Part I: Theoretical considerations. Geophysics , 62(5): 1586–1595. | |
| Berkhout A J, Verschuur D J. 2003. Transformation of multiples into primary reflections.//73rd Annual International Meeting, SEG, Expanded Abstracts, 1925-1928. | |
| Berkhout A J, Verschuur D J. 2004. Imaging multiple reflections, the concept.//74th Annual International Meeting, SEG, Expanded Abstracts, 1273-1276. | |
| Berkhout A J, Verschuur D J. 2006. Imaging of multiple reflections. Geophysics , 71(4): S1209–S1220. | |
| Berkhout A J, Verschuur D J. 2011. Full wavefield migration, utilizing surface and internal multiple scattering.//81th Annual International Meeting, SEG, Expanded Abstracts, 3212-3216. | |
| Brown M P, Guitton A. 2005. Least-squares joint imaging of multiples and primaries. Geophysics , 70(5): S79–S89. | |
| de Bruin C G M, Wapenaar C P A, Berkhout A J. 1990. Angle-dependent reflectivity by means of prestack migration. Geophysics , 55(9): 1223–1234. | |
| Engl H W, Hanke M, Neubauer A. Regularization of Inverse Problems. Boston: Kluwer Academic Publishers, 1996 . | |
| Fomel S. 2009. Adaptive multiple subtraction using regularized nonstationary regression. Geophysics , 74(1): V25–V33. | |
| Guitton A. 2002. Shot-profile migration of multiple reflections.//72nd Annual International Meeting, SEG, Expanded Abstracts, 1296-1299. | |
| Guo S J, Li Z C, Tong Z Q, et al. 2011. Joint imaging of primaries and surface-related multiples based on generalized shot-profile migration. Chinese J. Geophys. (in Chinese) , 54(4): 1098–1105. doi: 10.3969/j.issn.0001-5733.2011.04.025. | |
| Guo S J, Li Z C, Tong Z Q, et al. 2012. Method and technique for imaging of surface-related multiples. Progress in Geophys. (in Chinese) , 27(6): 2570–2576. doi: 10.6038/j.issn.1004-2903.2012.06.034. | |
| Hampson D. 1986. Inverse velocity stacking for multiple elimination.//56th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 422-424. | |
| He R Q, Schuster G. 2003. Least-squares migration of both primaries and multiples.//73rd Annual International Meeting, SEG, Expanded Abstracts, 1035-1038. | |
| Li Z C, Zhang J H. Seismic Data Processing Methods (in Chinese). (in Chinese) Dongying: China University of Petroleum Press, 2004 . | |
| Liu Y K, Chang X, Jin D G, et al. 2011. Reverse time migration of multiples for subsalt imaging. Geophysics , 76(5): WB209–WB216. | |
| Ryu J V. 1982. Decomposition (DECOM) approach applied to wave field analysis with seismic reflection records. Geophysics , 47(6): 869–883. | |
| Shan G J. 2003. Source-receiver migration of multiple reflections.//73rd Annual International Meeting, SEG, Expanded Abstracts, 1008-1011. | |
| Shan G J, Guitton A. 2004. Migration of surface-related multiples: tests on the Sigsbee2B dataset.//74th Annual International Meeting, SEG, Expanded Abstracts, 1285-1288. | |
| Taner M T, Koehler F. 1969. Velocity spectra-digital computer derivation applications of velocity functions. Geophysics , 34(6): 859–881. | |
| Tikhonov A N. 1963. Solution of incorrectly formulated problems and the regularization method. Soviet Mathematics-Doklady , 4(4): 1035–1038. | |
| Verschuur D J, Berkhout A J. 1997. Estimation of multiple scattering by iterative inversion, Part II: Practical aspects and examples. Geophysics , 62(5): 1596–1611. | |
| Weglein A B, Gasparotto F A, Carvalho P M, et al. 1997. An inverse-scattering series method for attenuating multiples in seismic reflection data. Geophysics , 62(6): 1975–1989. | |
| Wiggins J W. 1988. Attenuation of complex water-bottom multiples by wave-equation-based predication and subtraction. Geophysics , 53(12): 1527–1539. | |
| 郭书娟, 李振春, 仝兆岐, 等. 2011. 基于广义的炮偏移方法实现地表多次波和一次波联合成像. 地球物理学报 , 54(4): 1098–1105. | |
| 郭书娟, 李振春, 仝兆岐, 等. 2012. 表层多次波成像方法技术研究. 地球物理学进展 , 27(6): 2570–2576. | |
| 李振春, 张军华. 地震数据处理方法. 东营: 中国石油大学出版社, 2004 . | |