当今的油气资源严重供不应求，这就需要在老油田、老区块进行由简单的构造勘探向寻找复杂的隐蔽岩性油气藏转变的同时，还需要向新的未勘探过地区迈进，尤其复杂地表地区.剧烈起伏地表给地震资料采集、处理带来了严重的影响，研究者发展了一系列方法来克服这些问题(李振春等，2007；叶月明等，2008；刘国峰等，2010；张东良等，2011；张浩和张剑锋，2012).目前，针对起伏地表的处理主要有两种策略：一是对表层波场进行校正(Berryhill，1979，1984；李录明和罗省贤，2001)；二是基于起伏地表进行深度偏移成像(王成祥等，2002；安琪等，2007；杨海生，2007；雷秀丽等，2012).但复杂条件下静校正量难以准确计算，且静校正无法彻底消除地表起伏对地震波场造成的扭曲，因此直接基于起伏地表的深度偏移成像方法逐渐成为研究的热点.

随着计算机技术的发展，逆时偏移成为地震偏移成像领域重要的手段.Hemon(1978)详细研究了利用FD(Finite Discrete)求解波动方程的方法；随后Baysal等(1983)用Pseudo-spectral法实现了逆时偏移；Whitemore(1983)也提出了逆时偏移的研究思路，展示了逆时偏移的优越性；和LoewenthalMufti等(1983)通过利用指数解析解对双程波动方程进行求解，将逆时偏移的思想引入到地震勘探中.近年来，逆时偏移主要向叠前方向发展.叠前RTM的研究首先从声波方程开始的.Lesage等(2008)将有限差分和伪谱法相结合求解双程声波方程；Jones等(2007)将逆时偏移用来对复杂构造进行成像，并证实了逆时偏移对复杂介质成像的优势；徐义(2008)提出格子法叠前逆时偏移，适用于起伏地表数据的叠前逆时深度偏移成像.但是随着研究的深入，越来越多的学者对弹性波逆时偏移开展了大量的研究工作.Chang和Mcmechan(1987)通过FD方法实现了弹性逆时偏移；尧德中(1994)研究了单程弹性波逆时偏移，并与相移偏移进行了分析对比；Huseyin(2008)对波场进行分离之后再进行逆时偏移，这样可以有效地消除传统逆时偏移方法的成像假象，但同时增加的计算量，在处理复杂介质时，波场分离可能会引入大量的噪声，从而降低成像质量.Zhang和Zhang(2009)采用一步延拓法实现逆时深度偏移，该方法利用平方根算子将双程波动方程转化成类似单程波动方程的一阶偏微分方程形式，提高了逆时偏移的计算效率；但对基于波动方程的起伏地表逆时偏移研究很少，而且基本以声波为主(Zhang et al.，2010；王延光和匡斌，2012；孙小东等，2012).

由于矩形网格剖分对处理起伏地表存在的天然劣势，国内外地球物理学者对不规则网格地震波模拟与成像做了大量的研究.Fornberg(1988)在新的坐标系下应用伪谱法模拟波场，在模拟过程中直接基于曲网格进行；(Tessmer et al.，1992;Tessmer and Kosloff，1994)，Heshtolm和Ruud(1994)，Nielsen等(1994)在曲坐标系矩形网格下应用弹性波波动方程，用一个辅助坐标系将不规则起伏区域变换为规则区域.赵景霞等(2003)利用块映射和超限插值技术将曲界面变换成曲坐标系下的水平界面，并在新坐标下利用伪谱法模拟波场.王祥春和刘学伟(2005)利用双平方根法将地表采集的波场在曲坐标系下延拓，去除了起伏地表带来的影响.褚春雷和王修田(2005)等将三角网格应用于地震正演模拟中，能适应地表起伏的不规则性，但生成网格需要大量的计算量.

综上所述，起伏地表逆时偏移基本停留在声波介质的研究上，而国内外学者暂无在辅助坐标系下进行弹性波逆时偏移的尝试.

本文在传统坐标变换方法的基础上提出了分层坐标变换方法，并将此方法应用到了弹性波正演模拟和叠前逆时偏移中.在编程实现算法的基础上，通过对简单起伏模型和典型中原起伏模型进行正演模拟和成像试算，来验证本方法的正确性和实用性，并通过两个极端起伏地层高程不准确的情况进行测试来验证本方法的适用性.

传统坐标变换法将描述物理域直角坐标系下(x，z)的曲网格映射为计算域曲坐标系下(ξ，η)的正交网格(图 1A).但此方法在将起伏地表变为水平地表的同时会破坏原有的地下构造，对模拟和成像的结果造成一定的影响.本文提出的分层坐标变换法在传统坐标变换法的基础上进行改进，分别将每一层在物理域离散为曲网格，并将每一个起伏界面映射为水平界面(图 1B).分层坐标变换法所采用的映射函数为：

图 1

Fig. 1

图 1 网格剖分示意图(A)传统坐标变换；(B)分层坐标变换；a1，b1：曲网格坐标系；a2，b2：矩形网格坐标系. Fig. 1 Diagram of discrete grids(A)Traditional mapping method;(B)Layered mapping method； a1，b1: Curved grid system; a2，b2: Rectangular grid system.

由方程(2)可知，增加[(z_i-1(ξ)-z_i(ξ)]的值或减小(η_i-1-η_i)的值使(η_i-1-η_i)＜min[(z_i-1(ξ)-z_i(ξ)]都可满足稳定性，为了保全原始物理模型，通常减小(η_i-1-η_i)而不是增大[(z_i-1(ξ)-z_i(ξ)].当i=1，η₁=0，z₁(ξ)=0时，(1)、(2)式变为传统坐标变换方法.

曲网格上的参数值由实际模型插值得到.坐标变换仅仅是垂向上的拉伸或者压缩，采样点在坐标变换前后不发生变化，换句话说，曲网格上的波场值与矩形网格上的波场值是一一对应的.

当曲网格被映射为矩形网格时，已知一阶速度应力弹性波方程为：

坐标变换法的优势在于处理起伏自由地表. 本 文采用Hestholm and Ruud(2000)推导出的坐标 变换边界条件(6)式，该方法处理地表灵活而且具有较好的稳定性.

用二阶交错网格有限差分求解边界条件(6)式.然后随深度的增加，每个网格点增加两阶，直到十阶.时间上用二阶差分方法.通过采用高阶差分法可以提高计算精度，减少数值频散，也可减少地表地形的不规则性对体波和面波的散射.

对(6)式进行离散：

逆时偏移的第二步是沿时间反向计算检波点的检波波场.坐标变换法同样可以应用于波场逆推的过程中.采用跟坐标变换法正演一致的波场传播过程(式3)，从最大时刻t_max开始进行波场逆推，得到地下各个时刻的逆时外推波场.

利用亥姆霍兹分解方法分离纵横波波场是利用纵波是无旋场，横波是无散场的原理，近似得到弹性 波波场的散度为纵波波场，旋度为横波波场(Dellinger and Etgen，1990；Sun，1999).该方法原理比较简单，计算成本低，且计算精度在误差允许的范围内，因此常被用作分离波场.在正向和反向波场计算过程中，采用亥姆霍兹分解方法对震源波场和检波波场进行分离，本文推导得到分层坐标变换法计算域曲坐标系下的分离后二维波场近似表达式：

本文中采用的成像条件为震源归一化互相关成像条件(Kaelin and Guitton，2006)，即用震源照明强度的平方对运用互相关成像条件得到的成像结果进行归一化，得到归一化成像结果.对波场分离过的P波和S波应用互相关成像条件以减少每一分量P波和S波波场耦合在一起而造成的波场串叠现象.入射波场是P波，震源波场中P波能量远大于S波，所以用震源波场中的S波归一化，效果不如用P波归一化更接近正确的反射振幅.由此可知，P波入射时，用震源中的P波波场归一化，S波入射时，用震源中的S波波场归一化更合理准确.震源归一化互相关成像条件的数学表达式为：

该成像条件也称为应用照明补偿的互相关成像条件，它等价于应用频域匹配滤波的反褶积成像条件.它有效地消除了震源强度对成像结果的影响，成像振幅更接近与真实的反射系数.另外，还可以有效地补偿由于地震波在地层中传播时会产生能量损失，有利于深部地层成像.

对于逆时偏移还存在一个非常重要的问题，就是如何处理成像过程中产生的低频噪声.本文中采用Laplacian滤波方法(杜启振等，2013)，该方法能 在滤去低频噪声的同时，保证有效信号基本不受影响.

使用坐标变换逆时偏移方法需要注意的是，在进行逆时偏移的时候，输入的初始速度场为坐标变换之后的速度场.在逆时偏移过程中有两种思路：(一)波场的正向延拓、逆时延拓、波场分离及互相关成像都是在曲坐标系下完成的，得到的成像结果为曲坐标下的像，需要进行反变换得到最终的成像结果；(二)在每一时刻将在计算域得到的正向延拓和逆时延拓的波场值反变换到物理域，然后在物理域做互相关，最终叠加得到起伏地表成像结果.

比较两种思路，第二种思路的计算量略大于第一种思路，但是当速度模型较小，网格间距较大时，在反变换成像结果过程中，因为分辨率和反变换过程中线性插值的误差，成像结果会出现阶梯形状的同相轴断裂，对成像结果造成影响.为了消除此现象，网格间距需要加密很多，但是因此会造成计算量和内存的数倍增加.因此本文采用第二种思路.起伏地表坐标变换法逆时偏移的流程图如图 2所示.

图 2

Fig. 2

图 2 坐标变换法逆时偏移流程 Fig. 2 The flow chart of mapping RTM

另外，在应用此方法之前，需要得到近地表附近的地层高程，这可以通过传统有限差分方法或者其他起伏地表的成像方法得到.然后将地表附近的起伏底层离散成曲网格，而深部构造采用矩形网格.起伏地表对地下深部构造的影响要远大于地表的影响，因此相比于深部构造，得到浅层构造位置形态要容易得多.同时本文通过对两种极端情况进行测试可知，即使地表附近起伏高程不准确，甚至只对地表区域进行坐标变换以消除起伏地表的影响而对地下区域不变换以避免将地下复杂构造破坏，所得的成 像结果远好于传统坐标变换法逆时偏移的成像结果.

首先对图 3所示的起伏地表层状模型进行试算.图 3a为物理域起伏地表模型速度场，速度场的大小为4350 m×2840 m，起伏地表以上的速度为0，从上至下各层的纵波速度分别为2000 m·s^-1、 3000、3500、4000 m·s^-1，纵横波速度比约为1.732.横向采样点n_x=870，纵向采样点数n_z=568，网格间距为5 m×5 m，震源地表附近中点激发，震源主频为25 Hz，检波器位于地表，道间距为10 m.时间采样步长为0.2 ms，总记录时间为1.5s.图 3b为采用本文提出的分层坐标变化法得到的坐标变换之后的速度场；图 3c为采用传统坐标变换法对原速度场的参数进行插值计算，得到坐标变换后的速度场，网格离散及坐标变换示意图如图 4所示.从图 4可以看出：传统坐标变换虽然可以将起伏地表映射为水平地表，但地下构造也被破坏，影响了网格离散的精度；而本文提出的分层坐标变换法通过将起伏地表和地层分别映射为水平界面，很好的适应了有限差分矩形网格剖分的特点.分别采用本文提出的分层坐标变换方法和传统坐标变换方法进行试算，得到两种方法不同时刻的波场快照如图 5所示，得到的炮记录如图 6所示.

图 3

Fig. 3

图 3 起伏地表模型图(a)曲网格坐标系；(b)分层矩形网格坐标系；(c)不分层矩形网格坐标系. Fig. 3 Irregular surface model(a)Curved grid coordinate system;(b)Layered rectangular grid coordinate system;(c)Unlayered rectangular grid coordinate system.

图 4

Fig. 4

图 4 离散网格示意图(A)分层坐标变换法；(B)传统坐标变换法；(a1，b1)坐标变换前；(a2，b2)坐标变换后. Fig. 4 Diagram of discrete grids(A)Layered mapping method;(B)Traditional mapping method;(a1，b1)Before coordinate transform;(a2，b2)After coordinate transform.

图 5

Fig. 5

图 5 不同时刻波场快照图(A)分层坐标变换；(B)传统坐标变换；a1，a3，b1，b3：水平分量；a2，a4，b2，b4：垂直分量; a1，a2，b1，b2: 600 ms; a3，a4，b3，b4: 800 ms. Fig. 5 Snapshot of different time(A)Layered mapping method;(B)Traditional mapping method; a1，a3，b1，b3: Horizontal component; a2，a4，b2，b4: Vertical component; a1，a2，b1，b2: 600 ms; a3，a4，b3，b4: 800 ms.

图 6

Fig. 6

图 6 弹性波炮记录(A)分层坐标变换法；(B)传统坐标变换法；a1，b1：水平分量；a2，b2：垂直分量. Fig. 6 Elastic shot records(A)Layered mapping method;(B)Traditional mapping method; a1，b1: Horizontal component; a2，b2: Vertical component.

对比图 5所示的波场快照及图 6所示的炮记录，无论是传统坐标变换方法还是本文提出的分层坐标变换方法都能较好地处理起伏地表，入射波、反射波、转换波、面波等都可以清楚地模拟出来，因此坐标变换法对于处理起伏地表具有更好的效果.但同时也发现，本文提出的分层坐标变换方法除了在处理起伏地表方面有较好的效果外，相比传统坐标变换法来说，地下起伏地层的模拟也更精确，反射同相轴清晰，没有传统矩形网格阶梯状“毛刺”造成的虚假反射.

分析图 6A所示的采用本文提出的分层坐标变换法得到的炮记录可以看出：顶边界采用自由边界时，虽然采用的是纵波震源，但是由于自由边界的存在，炮记录中出现了直达S波(震源在地表发生反射产生的转换波)和还有面波，炮记录中面波与直达S波几乎重叠，使得直达S波同相轴较粗，不易看到.除了这两种波，地下反射层产生的反射波传播到自由地表发生反射，产生地表反射波和地表转换波.

为了更好的对比分析本文提出的分层坐标变换方法的优越性，本文对分层坐标变换法和传统坐标变换法得到的炮记录分别抽取了第120道波形，得到波形图 7A，图中各种波形如下：①直达P波；②面波和直达S波；③第一层的反射P波；④第二层反射P波；⑤第三层反射P波；⑥第一层转换PS波；⑦第二层转换PS波.放大虚框所示的区域得图 7B，本文提出的分层坐标变换法没有矩形网格剖分引起的虚假反射.通过波形图的对比可以看出本文提出的分层坐标变换方法对起伏地表和地下起伏地层的模拟都有很高的精度.

图 7

Fig. 7

图 7 (A)从图 6所示的炮记录中抽取的第120道波形图；(B)局部放大图a1，a3，b1，b3：水平分量；a2，a4，b2，b4：垂直分量；a1，a2，b1，b2：传统坐标变换法；a3，a4，b3，b4：本文提出的分层坐标变换法. Fig. 7 (A)The 120 trace waveform extracting from the shot records of Fig. 6;(B)Partial enlarged figure a1，a3，b1，b3: horizontal component; a2，a4，b2，b4: vertical component; a1，a2，b1，b2: traditional mapping method; a3，a4，b3，b4:layered mapping method proposed in the paper.

图 8为采用本文提出的分层坐标变换起伏自由 地表边界条件得到的波场分离快照图.从图中我们可以看到，纵横波被完全清晰地分离开来，包括入射波、反射波以及地表反射波和转换波.

图 8

Fig. 8

图 8 分层坐标变换法波场分离快照(A)600 ms;(B)800 ms；a1，b1：P波；a2，b2：S波. Fig. 8 Wave-field separation snapshot using layered mapping method(A)600 ms;(B)800 ms; a1，b1: P-wave; a2，b2: S-wave.英文标题

下面对该起伏层状模型进行多炮弹性波逆时偏移试算.使用100炮主频为25 Hz的雷克子波震源在地表附近放炮，中间激发两端接收，炮间距20 m，每炮470道接收，道间距5m.对多炮记录去除面波并进行表层多次波压制，将此多炮记录按照分层坐标变换技术进行逆时延拓及波场分离，并与正向延拓的分离波场进行互相关得到的PP成像和PS成像结果如图 9所示.图中PP成像和PS成像的无论地表附近还是深部地层都成像位置准确，同相轴清晰，且信噪比高.说明本文提出的分层坐标变换方法是准确有效的.从图中可以看出PP成像能量强于PS成像，但由于反射P波比反射S波的反射角大，因此检波器接收到的反射S的范围更广，PS成像范围更大.

图 9

Fig. 9

图 9 多炮成像结果(a)PP成像；(b)PS成像. Fig. 9 The image results using multi-shots(a)PP image;(b)PS image.

中原模型(图 10)起伏高程最大高差为260 m.横向网格点数n_x=1740，纵向网格点数n_z=1000.时间采样间隔为0.5 ms，记录时间为3s.采用纵波震源，震源主频为25 Hz，中点激发两端接收，每炮135道接收，检波器位于起伏地表处，道间距为16 m，炮间距为48 m，共100炮.图 11a为第10炮所覆盖的局部速度场(红色虚线区域).图 12为地表附近四层的起伏高程，为了避免地下复杂构造的原始形态遭到破坏，本文对地表附近地层(蓝色虚线区域)进行曲网格剖分，对地下深部区域进行矩形网格剖分，并进行分层坐标变换.图 11b为分层坐标变换后的局部速度场，也是分层坐标变换法逆时偏移时输入的速度场.与传统坐标变换后的局部速度场(图 11c)相比，分层坐标变换后的速度场很好的保留了地下深部构造的原始形态，而传统坐标变换后地下复杂构造变换为更复杂的构造，将会对成像造成很多未知的结果，同时网格剖分过程中会出现很多新的绕射点，严重影响成像结果.图 13为分别采用本文提出的分层坐标变换法和传统坐标变换法得到的第10炮记录，从图中可以看出，采用本文提出的分层坐标变换方法反射同相轴清晰，几乎没有虚假反射，而传统坐标变换炮记录虚假反射和绕射较多，且同相轴出现不连续.将多炮记录应用零延迟互相关成像进行成像，并对该成像结果进行Laplacian滤波，得到PP成像和PS的成像结果(图 14).

图 10

Fig. 10

图 10 中原起伏模型 Fig. 10 Zhongyuan topography model

图 11

Fig. 11

图 11 (a)局部速度场；(b)分层坐标变换后；(c)传统坐标变换后 Fig. 11 (a)Local velocity;(b)Layered mapping;(c)Traditional mapping

图 12

Fig. 12

图 12 地表附近四层界面的高程.(a)全局高程；(b)局部高程. Fig. 12 The elevation of four interface near surface(a)Global elevation;(b)Local elevation.

图 13

Fig. 13

图 13 第10炮单炮记录(A)分层坐标变换法；(B)传统坐标变换法；a1，b1：水平分量；a2，b2：垂直分量. Fig. 13 The 10th shot record(A)Layered mapping method;(B)Traditional mapping method; a1，b1: Horizontal component; a2，b2: Vertical component.

图 14

Fig. 14

图 14 中原模型分层坐标变换法成像结果(a)PP成像；(b)PS成像. Fig. 14 The image result of Zhongyuan model using layered mapping method(a)PP image;(b)PS image.

从图中可以看到，Laplacian滤波后的成像结果低频噪声消除的较为彻底，PS成像的能量弱于PP成像的能量，特别是深层位置，但PS成像范围要比PP成像范围广.无论是PP成像还是PS成像，深层反射界面都得到很好的成像，而且起伏地表附近同相轴形态清晰，具有较高的信噪比，震源的干扰效应正好能够体现起伏地表的形态.

作为对比，本文采用传统坐标变换逆时偏移方法得到的成像结果如图 15所示.从成像结果对比中可以看出，采用传统坐标变换逆时偏移得到的成像结果在浅层同相轴出现轻微错动，构造形态出现误差，在深层成像能量较弱，而且信噪比较低.而采用本文提出的分层坐标变换逆时偏移方法的成像结果同相轴清晰且连续，深部成像能量高，具有很高的信噪比.因此 本方法的成像结果无论从成像精度还是信噪比方面都更高.

图 15

Fig. 15

图 15 中原模型传统坐标变换法成像结果(a)PP成像；(b)PS成像. Fig. 15 The image result of Zhongyuan model using traditional mapping method(a)PP image;(b)PS image.

因为在实际情况中，地表附近的起伏高程可能也很难准确获得，因此我们采用两种极端情况来测试分层坐标变换方法的适应性.第一种情况，如图 16a所示，地下三层均采用平层进行分层坐标变换，得到的分层坐标变换后的局部速度场如图 17a；第二种情况，如图 16b所示，只将该模型分成地表附近和地下深层两层进行分层坐标变换，对地表附近进行曲网格剖分，对地下深层进行矩形网格剖分，得到的分层坐标变换后的局部速度场如图 17b.通过对比17a、b与11c可以看出，即使采用两种极端的情况进行分层坐标变换速度场原始构造形态的破坏程度远小于传统坐标变换.两种极端情况分层坐标变换逆时偏移得到的成像结果如图 18所示，这里只展示PP成像.通过成像结果可以看出，采用分层坐标变换法的成像结果明显好于传统坐标变换方法.因此，应用分层坐标变换方法时，在未得到地表附近的地下高程时，可以对地表附近一定区域采用曲网格剖分，对地下深部采用矩形网格剖分.

图 16

Fig. 16

图 16 两种极端情况下分层坐标变换法使用的地层高程(a)第一种情况；(b)第二种情况. Fig. 16 The interface elevation of two extreme cases using layered mapping method(a)The first case;(b)The second case.

图 17

Fig. 17

图 17 两种极端情况下坐标变换后的局部速度场(a)第一种情况；(b)第二种情况. Fig. 17 The velocity after coordinate transformation in two extreme cases(a)The first case;(b)The second case.

图 18

Fig. 18

图 18 中原模型两种极端情况下分层坐标变换法成像结果(a)第一种情况；(b)第二种情况. Fig. 18 英文标题

本文改进了传统坐标变换方法，实现了分层坐标变换起伏地表正演模拟，将起伏地表的映射推广到了地下起伏地层构造.同时，实现了基于变坐标系思想的起伏地表弹性波逆时偏移成像方法，并将本文中实现的分层坐标变换方法引入到了弹性波逆时偏移中，在实现算法的基础上，对简单起伏地表模型、典型起伏中原模型进行了正演模拟和成像试算，得到如下几点认识：

1)分层坐标变换方法不仅可以将起伏地表映射为水平地表，而且地下起伏地层也可以映射为水平地层，这样很好的适应了有限差分方法矩形网格剖分的特点，既可以保留传统有限差分算法的优点，又可以弥补传统算法对处理自由起伏地表的缺点，具有很高的模拟精度.

2)本文将坐标变换方法应用于起伏自由地表弹性波逆时偏移中，模型试算的结果证明了该方法无论在地表附近还是深层都具有较好的成像效果，同时震源的干扰效应能准确的体现起伏地表的形态.

3)在起伏地表的正演模拟和成像过程中，采用分层坐标变换方法可以仅对地表附近的起伏地层应用分层坐标变换，对深部区域应用传统矩形网格剖分，这样既可以较好的处理起伏地表对模拟和成像效果的影响，也可以完整的保持地下构造形态.

4)当地表附近地下起伏地层的高程无法准确获得时，使用近似起伏高程进行分层坐标变换法逆时偏移时，其成像精度较传统坐标变换法也有明显的提高.

5)应用变坐标系的方法，在辅助坐标系下进行逆时偏移的思想在未来处理起伏地表复杂构造具有广阔的前景；同时分层坐标变换的方法也可以直接应用到最小二乘偏移、全波形反演等方法中.

虽然本文提出的分层坐标变换正演模拟和逆时偏移方法在处理起伏地表时具有很大的优势，但此方法还不是特别完善，还存在很多需要完善的地方： 1)在应用坐标变换法进行正演模拟和逆时偏移时，要求地表起伏不能太剧烈，地表高程函数单调光滑，所以在进行剧烈起伏地表模型的逆时偏移时会出现不稳定现象； 2)而分层坐标变换方法对适用条件的要求更苛刻，它要求各个分层映射的函数都单调光滑，因此，在一些剧烈起伏或者断层构造区域，要进行部分分层坐标变换(如3.2)；3)传统弹性波逆时偏移存在计算量大，占用内存多的缺点，在分层坐标变换法弹性波逆时偏移中，这两个缺点依旧存在，在后续研究中，将在改进算法的同时，将GPU算法和双变算法(张慧等，2011)应用于分层坐标变换法弹性波逆时偏移中，提高起伏地表逆时偏移成像效率.