2019, Vol. 62
常规海洋反射地震理论假设海面是没有海浪的水平镜面.然而，海浪的存在将导致海面上下起伏，形成具有随机不规则起伏特性的自由表面.在理论上，不规则起伏海面将对入射地震波形成复杂的散射，进而对地下靶区的二次地震照明产生不同程度的影响，使反射地震观测数据的运动学和动力学特征相对于水平海面而发生变化.随着时移地震的发展和对勘探精度要求的提高，起伏海面对反射地震成像结果的影响变得不可忽略.事实上，Laws和Kragh(2000)对观测得到的时移地震成像结果进行分析，认为起伏海面对时移地震成像结果引起的误差不可忽略；Laws和Kragh(2010)系统地分析了起伏海面对地震数据的影响；指出即使在适度平静的海面环境下，起伏海面对反射地震数据某些频段上的振幅、相位的影响依然较大；Orji等(2012)利用一种基于克希霍夫-亥姆霍兹数值模拟方法计算了海面的散射波场；Egorov等(2015)采用克希霍夫近似法进行起伏海面下的数值模拟，分析了震源较深的情况下起伏海面对地震成像的影响.
起伏海面条件下的偏移成像、数值模拟问题与起伏地表下的有较大差异.首先，地表的起伏往往是大尺度的，可采用目前最为有效的贴体网格剖分格式进行数值求解(例如：孙建国和蒋丽丽，2009；蒋丽丽，2010；刘志强，2017).而海面的小尺度随机起伏特征会导致贴体网格的生成结果不收敛.针对这个问题，本文采用能够处理小尺度不规则起伏界面的不等距有限差分法.其次，在成像过程中，起伏地表条件下的地震数据来自于与地形相关的不同高程，需要将不同高程地震数据通过不同的手段进行“平化”，以校正起伏地表对成像结果的影响(王延光和匡斌，2012；刘定进等，2016；侯爵等，2018).不同的是，起伏海面下的地震数据均来自于水下同一深度，海面影响主要体现在鬼波和自由表面散射上.需要“平化”的仅仅是地震数据中由海面产生的散射，而非起伏地表下不同高程所引起的到时误差.最后，由于海面形态的随机不确定性，导致了常规的“平化”的方法无法对起伏海面的影响进行有效消除.这意味着，研究起伏海面对成像结果的影响具有重要的理论和实际意义.
本文将重点研究随机起伏海面对波场照明和偏移成像结果的影响.利用海浪谱建立起伏海面模型，通过对波动方程的不等距差分离散，从地震波照明的角度分析了起伏海面对偏移成像结果的影响.在具体的技术路线上，本文首先利用不等距有限差分法生成含有起伏海面影响的反射地震数据，然后对所形成的数据进行反向外推和成像，最后对成像结果进行分析和对比.
1 基本原理为了研究起伏海面对地震波场传播的影响，国内外许多学者普遍使海面“冻结”到某一个时刻上，即“冻结”海面模型.综合前人的工作，本文选用能够处理小尺度起伏问题的不等距有限差分法进行起伏海面条件下波动方程的求解.
图 1a是不等距差分格式示意图.将网格角点上的点称为正则内点，不与角点重合的点称为非正则内点.其中：h为网格间距；θ为不同方向上包含不规则边界高程信息的参数，即：真实高程与网格间距的比值；
|
图 1 不等距差分格式示意图和波场切片 Fig. 1 The diagram of the unequal distance difference scheme and wave field slice |
对于声波方程：
|
(1) |
不规则界面处x方向波场的二阶导数为：
|
(2) |
则声波方程的不等距离散格式(时间二阶空间二阶)为
|
(3) |
在复杂边界处求解波动方程的过程中，为满足空气-海水界面的波场(u=0)的边界条件，通常假设非正则内点p1、p2、p3、p4和界面上方的点的波场值不为0，这样的点称为鬼点(Zhang et al., 2013；孙章庆等, 2012a, b).鬼点的值与空间中波场值无关，不改变波场值的空间分布，仅仅参与波场导数的计算.
图 1b为由该方法求解的某一时刻起伏海面下的地震波波场快照.可见，由于海面的起伏，二次反射的形态受到一定的影响，在波尾产生散射和畸变现象.
2 起伏海面条件下的照明分析起伏海面使地震波形成复杂的散射进而影响靶区波场照明，下面将以上文中不等距差分法为基础对波场照明进行分析、对比.
图 2a为平静海面下的二次反射传播方向示意图；图 2b为起伏海面下发散的二次反射传播方向示意图；图 2c为起伏海面下聚焦的二次反射传播方向示意图；图 2d为起伏海面下向一个方向偏折的二次反射传播方向示意图.
|
图 2 起伏海面对地下照明的影响示意图 Fig. 2 The diagram of the influence of the undulating sea surface on the underground illumination |
地下照明代表着波场的空间分布，故在某一点(i, j)地震波场地下照明可以表示为：
|
(4) |
其中：q(i, j)表示空间中某一位置的波场照明强度；u(i, j, t)为某一时刻t空间中(i, j)处的波场值；T为总的记录时间.
2.2 均匀海水模型地下照明图 3a为水平海面下均匀海水地下照明结果；图 3b、3c、3d分别为不同起伏海面下均匀海水地下照明与水平海面下照明结果差异；其中震源在水下10 m、横向4 km处、主频20 Hz.可知在不同形态的海面下，地震波场的地下照明分布是不同的；图 3b中震源上方的海面处于近似波谷的状态，这时地震波场的地下照明分布在一定范围内发散(照明形如图 2b)，震源两侧的照明变好，震源下方照明变差；图 3c中震源上方的海面处于近似波峰的状态，这时地震波场的地下照明分布在一定范围内聚焦(照明形如图 2c)，震源两侧的照明变差，震源下方的照明变好；图 3d中震源上方的局部海面处于向某一方向倾斜的状态，这时地震波场的地下照明分布也变得倾斜，在震源右侧范围内照明变好，左侧变差(照明形如图 2d).
|
图 3 起伏海面对均匀介质地下照明的影响 Fig. 3 Influence of undulating sea surface on underground illumination of homogeneous medium |
图 3d、4a、4b分别为震源深度(10 m、20 m、30 m)同一起伏海面下均匀海水地下照明结果与水平海面地下照明结果的差；可见，震源在不同的深度时，地震波场的地下照明分布是不同的，且随着震源深度的增加，震源上方更大范围的起伏海面使地震波形成复杂的散射，导致起伏海面对波场照明的影响更复杂.
|
图 4 不同震源深度下起伏海面对均匀介质地下照明的影响 Fig. 4 Influence of undulating sea surface on underground illumination of homogeneous medium by the source at different depths |
图 5a为水平海面条件下的Marmousi模型(图 8)的地下照明结果；图 5b、5c、5d分别为不同起伏海面Marmousi模型地下照明与水平海面Marmousi模型地下照明结果的差异；受起伏海面影响，图 5b中震源两侧一定范围内照明变好，正下方照明变差；图 5c中震源两侧一定范围内照明变差，正下方照明变好；图 5d中震源左侧一定范围内照明变差，右侧照明变好.
|
图 5 地震波照明分析 Fig. 5 Seismic wave illumination analysis |
|
图 8 Marmousi速度模型 Fig. 8 The Marmousi velocity model |
综合上述所有的照明分析可知：起伏海面会改变地震波场地下照明，且对于不同的地质模型，同一起伏海面对波场照明的影响具有一定的相似性.
3 起伏海面下RTM成像起伏海面使波场靶区照明变得不均匀，进而影响观测数据的照明分布.这种照明不均在地震数据上的表现为振幅和相位的畸变.那么，将畸变的数据进行波场反推和成像，势必会产生一定的误差与噪声.为了研究起伏海面对成像结果的影响，本文以RTM方法为例进行成像、对比、分析.由于RTM方法为常规方法，故不在此详细论述，详见罗宾(2012).
具体流程上，首先以不同海面为变量，通过正演分别获得平静海面和起伏海面下的反射地震数据.然后分别对地震数据进行反推、成像，并对成像结果进行分析、对比.
3.1 起伏海面条件下简单模型的RTM成像为了研究起伏海面对空间横向不同位置成像结果的影响，建立图 6a所示的速度模型，该模型横向401个采样点，纵向301个采样点，网格间距5 m.
|
图 6 单层RTM成像 Fig. 6 Single layer model RTM imaging |
图 6b、6c分别为平静海面和起伏海面下该模型的单炮RTM成像结果.可见，由于起伏海面的影响，成像结果的噪声增多且成像界面发生弯曲畸变；对于噪声增多的现象，是由于用于成像的正向波场与反推波场中起伏海面散射信息相关的结果；对于成像界面的畸变和弯曲的现象，是由于用于成像的不含海面散射信息的正向波场与含有海面散射信息的反推波场相关的结果.图 6d为起伏海面下该模型的多炮RTM成像结果.随着照明条件的变好，成像界面的弯曲现象被“平均”，且变得不明显；而由起伏海面产生的成像噪声则在空间中多炮累加在一起，依旧存在.
为了研究起伏海面对模型不同深度的成像结果的影响，我们建立图 7a所示的水平多层模型.该模型横向401个采样点，纵向301个采样点，网格间距5 m，共五层，随着模型深度的增大，地震波的速度增加.
|
图 7 多层RTM成像 Fig. 7 Multi-layer model RTM imaging |
图 7b为起伏海面下该模型的单炮RTM成像结果.在较浅的深度上，成像结果噪声较多、成像界面弯曲和畸变现象较严重；随着深度增加和波场畸变扩散，成像结果噪声减少、成像界面弯曲和畸变现象减弱；对于深度较大的地层，成像噪声、成像界面的弯曲和畸变现象几乎不存在.
3.2 起伏海面条件下Marmousi模型的RTM成像上述主要通过简单层状模型的RTM成像结果来分析起伏海面对地震成像结果的影响.对于复杂模型，我们以Marmousi模型为例进行试算、分析.
图 8为Marmousi模型，该模型横向1601个采样点，纵向301个采样点，网格间距5 m.观测排列共有53炮，震源间隔150 m，采用时间二阶空间四阶差分算子，波场切片存储间隔为5 ms，采用完全匹配层(Perfectly Matched Layer, PML)吸收边界条件；此外，对所有的成像结果都进行Laplace滤波去除低频噪声.
图 9a为平静海面下Marmousi模型的RTM成像结果，图 9b、9c分别为不同起伏海面下Marmousi模型的RTM成像结果.经对比分析，可见起伏海面条件下成像结果整体上更“粗糙”，噪声较多.图 9d、9e、9f分别为平静海面下和不同起伏海面下模型A区域的成像结果放大图.可见，在起伏海面的成像结果中，含有主要沿着射线路径近似分布的波纹状成像噪声(图 9e、9f圈中).这种成像噪声主要是由于反射地震数据中包含的起伏海面散射信息经过反推和相关成像所引起的；此外，起伏海面条件下震源噪声无规律，这从侧面证实了起伏海面对波场地下照明的影响.与平静海面成像结果相比，在图 9e中，箭头处的成像界面向上发生畸变，而在图 9f中，箭头处的成像界面向下发生畸变；可知对于不同的起伏海面，成像界面可能发生不同方向的弯曲畸变.图 9g、9h、9i分别为平静海面和不同起伏海面下B区域的成像结果放大图.在平静海面成像结果图 9g圈中，可轻易分辨出该处的多个地层；而在起伏海面成像结果图 9h、9i圈中，多个地层混为一层，这是由于起伏海面散射引起的噪声和成像界面的畸变、弯曲经过多炮叠加后的结果.可知，起伏海面会降低成像结果的分辨率.此外在起伏海面成像结果(图 9h、9i)的箭头处，成像界面弯曲的现象依旧存在，且对于不同起伏海面，弯曲形态存在一定的差异.综上，对于复杂模型下的时移地震成像，起伏海面引起的散射会使成像界面发生畸变、使复杂构造处的成像分辨率降低.
|
图 9 不同海面下Marmousi模型的逆时偏移成像结果 Fig. 9 Inverse time migration imaging results of Marmousi models under different sea surface |
本文利用海浪谱建立起伏海面模型，对起伏海面条件下的地震波进行了数值模拟和照明分析，并在此基础上对起伏海面下的简单和复杂模型进行了成像结果的对比.得到下列结论和认识：
(1) 起伏海面对地震波场的地下照明分布产生一定的影响，这种影响主要与震源附近的海面形态有关，且不同形态的起伏海面对地震波场地下照明的影响特征不同.
(2) 相对于平静海面，起伏海面下单炮反射地震成像结果噪声增多，且成像界面发生畸变和弯曲；对于多炮反射地震成像，随着照明条件变好，成像界面的弯曲和畸变效应在一定程度上被平均，变得不明显，但由起伏海面引起的成像噪声依旧存在.
(3) 纵向上，起伏海面对于浅层的反射地震成像结果影响较大；横向上，起伏海面对于震源正下方或附近的反射地震成像结果影响较大.
(4) 对于复杂的地质模型，起伏海面使成像噪声增多且成像界面发生畸变和弯曲.一定程度上，降低了成像结果的分辨率，影响地下地质构造的精细解释.
综上所述，对于精度要求较高的时移地震，在某些情况下，为得到更满意的地震成像结果，起伏海面的影响不可忽略, 需要对地震数据或成像剖面进行多道信号特征平均的整形滤波，校正由起伏海面产生的振幅和相位畸变.
附录A 冻结的起伏海面建模海浪的存在使海面的形态呈现随机时变的特点.在20世纪40年代，海洋学界把起伏的海面形态看成由不同频率不同相位的波动所叠加的结果.在物理上，S(ω)可以看成频率为ω的海浪的能量，由S(ω)和ω构成的用于描述海水能量分布的谱即为海浪谱(邹建武等，2010).
一个区域的海浪谱通常是在海浪的不同成长时期和不同海洋环境(如风速、重力、水温等)下，经过多次观测所拟合的经验性结果.至今应用最多海浪谱的仍是PM(Pierson-Moscowitz)谱，其表达式为：
|
(A1) |
其中：ω为海面的起伏频率，u为海面上19.5 m处的风速，g为当地的重力加速度，α, β为常数，α=0.0081，β=0.74.
PM谱是由Pierson和Moscowitz在20世纪50年代对北大西洋充分成长型海浪进行观测所提出的海浪谱(Pierson and Neumann, 1955；Pierson and Moskowitz, 1964).由于PM谱的观测选取了海面上空19.5 m处的风速数据，所以PM谱中的参数u代表海面上19.5 m处的风速；对于不同的海浪谱，u代表不同高度的风速，且不同高度的风速之间满足一定的换算关系.对于重力加速的g，本文选取常用的9.8 m·s-2，除极端的纬度较高地区外，g的值影响不大.
|
附图 1 基于PM谱的起伏海面建模 Appendix Fig. 1 Modeling of undulating sea surface based on PM spectrum |
Egorov A, Glubokovskikh S, Bona A, et al. 2015. Influence of rough sea surface on sea surface reflections: deep towed high-resolution marine seismic case study.//85th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 3661-3665.
|
Hou J, Liu Y S, Lan H Q, et al. 2018. Elastic reverse time migration using a topography flattening scheme. Chinese Journal of Geophysics (in Chinese), 61(4): 1434-1446. DOI:10.6038/cjg2018L0624 |
Jiang L L. 2010. Body-fitted grid generation for the Geological conditions[Ph. D. thesis] (in Chinese). Changchun: Jilin University.
|
Laws R, Kragh E. 2000. Time-lapse seismic and the rough-sea wavelet.//70th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 19(1): 1603.
|
Laws R, Kragh E. 2010. Rough seas and time-lapse seismic. Geophysical Prospecting, 50(2): 195-208. |
Liu D J, Jiang B, Li B, et al. 2016. Rugged-topography reverse time migration in complex piedmont zone. Oil Geophysical Prospecting (in Chinese), 51(2): 315-324. |
Liu Z Q. 2017. Study on seismic wave simulation based on orthogonal curvilinear mesh under the complex geological conditions[Ph. D. thesis] (in Chinese). Changchun: Jilin University.
|
Orji O C, Söllner W, Gelius L J. 2012. Effects of time-varying sea surface in marine seismic data. Geophysics, 77(3): 33-43. |
Pierson W J, Neumann G, James R W. 1955. Practical methods for observing and forecasting ocean waves by means of wave spectra and statistics. Washington DC: US Naval Oceanographic Office.
|
Pierson W J Jr, Moskowitz L. 1964. A proposed spectral form for fully developed wind seas based on the similarity theory of S. A. Kitaigorodskii. Journal of Geophysical Research, 69(24): 5181-5190. DOI:10.1029/JZ069i024p05181 |
Sun J G, Jiang L L. 2009. Orthogonal curvilinear grid generation technique used for numeric simulation of geophysical fields in undulating surface condition. Oil Geophysical Prospecting (in Chinese), 44(4): 494-500. |
Sun Z Q, Sun J G, Han F X. 2012a. Traveltime computation using the upwind finite difference method with nonuniform grid spacing in a 3D undulating surface condition. Chinese Journal of Geophysics (in Chinese), 55(7): 2441-2449. DOI:10.6038/j.issn.0001-5733.2012.07.028 |
Sun Z Q, Sun J G, Han F X. 2012b. The comparison of three schemes for computing seismic wave traveltimes in complex topographical conditions. Chinese Journal of Geophysics (in Chinese), 55(2): 560-568. DOI:10.6038/j.issn.0001-5733.2012.02.018 |
Wang Y G, Kuang B. 2012. Prestack reverse-time depth migration on rugged topography and parallel computation realization. Oil Geophysical Prospecting (in Chinese), 47(2): 266-271. |
Zhang D L, Zhan G, Schuster G T. 2013. Multi-source least-squares reverse time migration with topography.//83rd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 3736-3740.
|
Zou J W, Zhu M B, Dong W. 2010. An overview of ocean wave modeling. Ship Electronic Engineering (in Chinese), 30(11): 10-14. |
侯爵, 刘有山, 兰海强, 等. 2018. 基于起伏地形平化策略的弹性波逆时偏移成像方法. 地球物理学报, 61(4): 1434-1446. DOI:10.6038/cjg2018L0624 |
蒋丽丽. 2010.面向地质条件的贴体网格生成技术[博士论文].长春: 吉林大学.
|
罗宾(Robein E.)[法]. 2012.地震资料叠前偏移成像——方法、原理和优缺点分析.石油工业出版社.
|
刘定进, 蒋波, 李博, 等. 2016. 起伏地表逆时偏移在复杂山前带地震成像中的应用. 石油地球物理勘探, 51(2): 315-324. |
刘志强. 2017.复杂地质条件下基于正交曲线网格的地震波数值模拟研究[博士论文].长春: 吉林大学
|
孙建国, 蒋丽丽. 2009. 用于起伏地表条件下地球物理场数值模拟的正交曲网格生成技术. 石油地球物理勘探, 44(4): 494-500. DOI:10.3321/j.issn:1000-7210.2009.04.020 |
孙章庆, 孙建国, 韩复兴. 2012a. 三维起伏地表条件下地震波走时计算的不等距迎风差分法. 地球物理学报, 55(7): 2441-2449. DOI:10.6038/j.issn.0001-5733.2012.07.028 |
孙章庆, 孙建国, 韩复兴. 2012b. 针对复杂地形的三种地震波走时算法及对比. 地球物理学报, 55(2): 560-568. DOI:10.6038/j.issn.0001-5733.2012.02.018 |
王延光, 匡斌. 2012. 起伏地表叠前逆时深度偏移与并行实现. 石油地球物理勘探, 47(2): 266-271. |
邹建武, 祝明波, 董巍. 2010. 海浪建模方法综述. 船舶电子工程, 30(11): 10-14. |