0 引言

中国中西部地区山前带断裂、高陡构造发育，地震反射数据信噪比低。相比于平缓地层，高陡构造区地震照明明显不足，特别是常规地震采集中各区域大多采用正交均匀炮检点布设方式，进一步凸显高陡构造区照明不足问题，严重影响成像品质。

针对高陡构造照明不足，影响叠加成像的难题，近年来许多学者已进行大量尝试，相关研究大致可归纳为以下四类。①采用更多类型的地震波场参与偏移成像，增加高陡构造区域的照明强度。如Guitton^[1]、Liu等^[2]、曲英铭等^[3]及秦宁等^[4]将多次波纳入有效信号范畴，研发了相应的多次波成像方法；Lu等^[5]根据散射波能量随散射角变化的特性，利用快速拉东变换成功分离散射波场，实现了基于散射波能量的地震成像。②提高偏移算子对大角度、高波数波场的成像精度，以充分利用来自高陡构造区域的信号，改善成像质量。Zhang等^[6]及Tang等^[7]通过优化波场延拓算子，有效提高傅里叶偏移方法对大角度波场的成像精度，取得很好的应用效果。③通过对成像条件进行优化，压制偏移噪声，突出大角度反射信号，改善成像质量。如逆时偏移(Reverse Time Migration，RTM)延拓得到的反射波场与正传波场互相关会产生较强低频噪声，这类噪声会严重干扰高陡界面的成像^[8-9]。Zhang等^[10]利用Laplace滤波压制此类低频噪声；王晓毅等^[11]通过Poynting矢量实现上下行波场的分离，并选择特定传播方向的波场做互相关，有效降低偏移噪声影响。这三类都是在偏移过程中改善高陡构造成像质量。④尝试从现场采集环节，增加目标位置照明强度，提高原始数据品质^[12-13]。该方案的典型实例包括基于点绕射成像理论的角度域照明分析^[14]和基于共聚焦偏移理论^[15-16]的照明分析^[17]。通过对采集观测系统的照明分析，不断调整观测系统设置和优选激发位置，提高地下地质目标的照明度。

上述四类方案在不同程度上改善了高陡构造的成像品质，但也存在一些不足，如“两宽一高”地震采集技术^[18-19]的广泛应用，海量数据偏移成像需成倍增加的计算及存储资源，前三类基于偏移过程实施的方案大多并未集成到商业软件平台，无法有效调配大规模计算资源应对海量数据的成像需求。对于第四种方案，高陡构造区域一般会伴生强起伏的地表条件^[20]，山崖峭壁区域无法正常开展地震数据采集。另外，照明分析依赖于速度模型，在勘探前期，很难获得复杂构造区高质量的速度模型，加之受采集成本限制，基于观测系统照明分析的方案最终实施效果常会大打折扣，难以满足高陡构造精确成像所需的照明强度。

与上述四种方案不同，本文建立了一种改善高陡构造成像质量的新方法。不同区域的高陡构造反射具有明显的方向性，在跨不同区域的长排列单炮地面记录上，不同炮检距(即入射角)、不同方位角段的反射信号具有明显的分带特征。首先根据长排列单炮偏移记录上不同炮检距、不同方位角段的成像特征进行分片，每个片受地震波传播角度控制，称为炮域炮检距矢量片(Shot domain offset vector tile，SOVT)，并将相关单炮的SOVT叠加获得SOVT剖面。随后基于共聚焦理论^[21]计算各个SOVT剖面的成像清晰度属性，并以此为依据对各个SOVT剖面进行优化叠加，突出有效信号、降低偏移噪声，在一定程度上校正因观测系统、频带限制等因素导致的成像误差，达到改善高陡界面成像质量的目的。本文方法主要对炮域深度偏移输出结果进行优化，可依托于商业软件中的深度偏移模块运行，对计算资源需求极低，适于工业化推广。

1 方法原理 1.1 炮域炮检距矢量片划分

在现场地震数据采集中，炮域数据是最常见的地震记录形式，记为U(X_s，X_r，T)，其中X_s为炮点位置，X_r为接收点位置，T为记录时间。炮域叠前深度偏移公式可表示为

$ I(X)=\sum\limits_{i=1}^n U_i\left(X_{\mathrm{s}}, X_{\mathrm{r}}, T\right) M_i $

(1)

式中：I(X)为X位置成像结果；U_i和M_i分别为偏移孔径内的炮域数据和偏移算子；n为偏移孔径内的炮数。采用式(1)输出偏移叠加剖面，地震传播方向信息(炮检距、方位角)将被叠加运算覆盖。

如引言所述，区域高陡构造反射具有明显方向性，反射信号常集中于特定炮检距和方位角范围内。常规深度偏移算法等权处理所有地震数据，不能突出信号优势方向特征，难以改善最终成像剖面的信噪比。基于此认识，在式(1)中加入地震激发传播方向项，在偏移孔径范围内，将单炮的深度偏移结果按特定炮检距和方位角输出并叠加，形成SOVT剖面，相关公式表示为

$ \operatorname{SOVT}(X)=\sum\limits_{i=1}^n K_i(\text { offset }, \theta) U_i\left(X_{\mathrm{s}}, X_{\mathrm{r}}, T\right) M_i $

(2)

式中：SOVT(X)为特定方向波场在X位置成像结果；offset为炮检距；θ为方位角；K_i为偏移孔径内单炮数据对应的地震激发传播方向项，筛选出单炮特定传播方向波场的偏移结果。

实际应用中，需根据地震—地质条件在偏移孔径内划分SOVT。当地层界面倾角与单炮地震波下传角度一致时，地表可接收到最大比例的反射信号，因此SOVT的炮检距划分应与地层倾角保持一致。

图 1左可围绕垂直入射、入射角α₁、入射角α₂划分SOVT的炮检距范围；同理，SOVT的方位角划分应与地层倾向保持一致，同时SOVT的数量应与地下构造的复杂程度相匹配。在图 1右的划分方案中，来自相对平缓界面的反射地震能量集中在SOVT1剖面，而来自高陡界面的地震能量分布在SOVT2~SOVT5中。SOVT方法的本质是基于地震波传播方向分别输出偏移成像结果，各类偏移算子(单程波、双程波等)皆可使用，可纳入当前主流偏移成像方法中(RTM、高斯束偏移等)。

1.2 基于共聚焦分析的优化叠加

常规深度偏移方法是对式(2)得到的所有SOVT剖面进行等权叠加。本文主旨是区别对待不同的SOVT剖面，压制噪声以突出高陡界面的弱有效信号，实现优化叠加，这需要定量分析各SOVT剖面的信号聚焦特征。本文利用共聚焦分析方法，分析各SOVT剖面的聚焦特性，提升含有效信号SOVT剖面的权重，实现高陡界面的高精度叠加成像。

Berkhout^[15]提出WRW模型，采用矩阵形式描述地震波从激发、反射到接收的全过程。地表Z₀处激发、接收的地震记录P(Z₀，Z₀)用WRW模型表示为

$ \begin{aligned} & \boldsymbol{P}\left(Z_0, Z_0\right)=\boldsymbol{D}\left(Z_0\right) \times \\ & \quad \sum\limits_{n=1}^N\left[\boldsymbol{W}\left(Z_0, Z_n\right) \boldsymbol{R}\left(Z_n, Z_n\right) \boldsymbol{W}\left(Z_n, Z_0\right)\right] \boldsymbol{S}\left(Z_0\right) \end{aligned} $

(3)

式中：D(Z₀)和S(Z₀)分别为接收和激发矩阵，包括接收、激发位置和相关特性(子波等)信息，地震波反射层由深度层Z_n表示；W(Z_n，Z₀)和W(Z₀，Z_n)分别代表目标深度Z_n和地表Z₀之间的下行和上行波传播矩阵；R(Z_n，Z_n)为目标深度Z_n处反射率矩阵。

当地震记录为单点绕射响应δP(Z₀，Z₀)时，可在接收端和激发端分别应用聚焦算子F(Z_n，Z₀)和F(Z_n，Z₀)消除地震采集因素及传播效应对地震波场的影响，得到绕射点的分辨率函数，相关公式为

$ \begin{aligned} {\rm{ \mathsf{ δ} }} \boldsymbol{P}_{j j}\left(Z_n, Z_n\right) & =\boldsymbol{F}\left(Z_n, Z_0\right) {\rm{ \mathsf{ δ} }} \boldsymbol{P}\left(Z_0, Z_0\right) \boldsymbol{F}\left(Z_0, Z_n\right) \\ & =\boldsymbol{B}_{D, j}\left(Z_n, Z_0\right) \boldsymbol{B}_{S, j}\left(Z_0, Z_n\right) \end{aligned} $

(4)

式中：下标j代表绕射点位于(X_j，Z_n)处；B_D，j(Z_n，Z₀)=F(Z_n，Z₀)D(Z₀)W_j(Z₀，Z_n)，代表接收端聚焦矩阵；B_S，j(Z₀，Z_n)=W_j(Z_n，Z₀)S(Z₀)×F(Z₀，Z_n)，代表激发端聚焦矩阵。B_D，j(Z_n，Z₀)与B_S，j(Z₀，Z_n)相乘即可得到绕射点的分辨率矩阵δP_jj(Z_n，Z_n)。该分辨率矩阵剔除了反射系数影响，可有效评估特定地震成像系统(观测系统、速度模型、频带宽度和偏移算子等)对任意点的成像能力，目前多用于指导地震观测系统的优化设计^[22]。

将式(2)输出的SOVT剖面进行加权融合是提高叠加质量的关键。本文采用共聚焦分析方法，评估各SOVT剖面局部位置的聚焦成像特征，以突出高质量成像单元。

首先在式(4)中加入传播角度项A_i(Z_n，Z_n)，计算不同入射角度波场的分辨率函数

$ {\rm{ \mathsf{ δ} }} \boldsymbol{P}_{j j, i}\left(Z_n, Z_n\right)=\boldsymbol{F}\left(Z_n, Z_0\right) \boldsymbol{D}\left(Z_0\right) \boldsymbol{W}_j\left(Z_0, Z_n\right) \boldsymbol{A}_i\left(Z_n, Z_n\right) \boldsymbol{W}_j\left(Z_n, Z_0\right) \boldsymbol{S}\left(Z_0\right) \boldsymbol{F}\left(Z_0, Z_n\right) $

(5)

式中下标i为SOVT编号。通过传播角度项A_i(Z_n，Z_n)可筛选特定传播方向的波场，其筛选方式应与SOVT输出方案一致。在实际计算中可基于SOVT划分方案的入射角、方位角信息，并结合埋深围绕成像点地表投影在地表划定各SOVT的对应区域完成波场筛选。通过式(5)可计算各SOVT剖面的成像分辨率函数，为后续优化叠加提供分析基础。

图 2为通过式(4)、式(5)得到的SOVT剖面成像分辨率函数(二维切片)，形态上类似于地震子波，包括主瓣和旁瓣两部分。一般认为主瓣聚焦越好、旁瓣扰动越小，则成像质量越好。理想状态下分辨率函数为脉冲，即只在聚焦点位置有数值，其余位置皆为零，代表成像系统能完全反映地下反射系数。通过分辨率函数可评估成像质量，但需同时兼顾主瓣和旁瓣的影响。

当两个相邻散射点的分辨率函数交点位于主瓣35%振幅处时，为地震水平分辨极限^[23]。因此，本文以大于35%振幅的主瓣代表成像有效能量，结合Schoenberger^[24]使用主瓣与旁瓣的峰值比评价成像分辨率的方法，求取有效成像能量在分辨率函数中的占比，得到成像清晰度(sharpness)

$ S=\frac{\sum \mathrm{AMP}_{35 \%}^2}{\sum\limits_{k=1}^M \sum\limits_{l=1}^N \mathrm{AMP}_{k l}^2} $

(6)

式中：AMP_35%代表大于35%的主瓣振幅；AMP_ij代表分析点周边的分辨率函数振幅，包括Inline和Crossline方向，M、N对应为此两方向采样点总数。成像清晰度表示偏移结果中有效能量占比，数值越大主瓣聚焦越好，旁瓣扰动越小，意味着更高成像质量，因此可通过成像清晰度属性评估SOVT剖面局部成像质量。

成像清晰度属性与地下构造和地表观测系统密切相关，对于变化剧烈区域，需加大分析密度，当样点数量足够控制空间变化特征时，可通过空间线性插值得到地下所有样点的成像清晰度属性，将其作为权叠加系数完成SOVT剖面的优化叠加。相关公式如下

$ I_{\mathrm{op}}\left(X_j\right)=\frac{1}{\sum\limits_{i=1}^m S_j^i} \sum\limits_{i=1}^m S_j^i \operatorname{SOVT}^i\left(X_j\right) $

(7)

式中：上标i为第i个SOVT剖面；下标j为X位置的第j个样点；S为成像清晰度；I_op(X_i)为经优化叠加后在X_i位置的成像结果。成像清晰度S的取值范围在0~1之间，当数值越大时，意味着成像质量越高，优化叠加的权重也越大。基于共聚焦分析的优化叠加运算除了可突出特定倾向、倾角界面的有效成像能量外，也可在一定程度上校正因观测系统分布不均、数据频带限制等因素导致的成像误差，改善高陡构造区成像质量。

基于共聚焦分析的SOVT剖面优化叠加方法实施过程细述如下：

步骤1：炮集常规深度偏移运算，得到各单炮偏移结果；

步骤2：根据区域构造特征的差异性及复杂程度，建立SOVT划分方案；

步骤3：集成各单炮SOVT，按照式(2)输出，得到多个SOVT剖面；

步骤4：基于共聚焦分析，按式(5)和式(6)求取不同SOVT剖面的单点成像清晰度属性，所求样点数量需足够控制地表和地下条件变化导致的成像清晰度属性变化，采用空间线性插值得到所有样点的成像清晰度属性；

步骤5：利用步骤4得到的成像清晰度属性，按式(7)完成所有SOVT剖面(步骤3结果)的加权叠加，得到最终优化成像剖面。

2 模型数据测试

为评价本文方法对高陡界面的成像效果，设计一个二维三层凹面速度模型(图 3)做正演测试。模型尺寸为2000 m×1500 m，纵横向采样间隔为10 m；设计的三层介质从上到下速度分别为2000、2500、4000 m/s，其中第一水平界面位于深度400 m处，第二界面位于深度800 m处，在第二界面中心位置存在一个直径为700 m的凹面，中心深度为200 m；所有激发点、接收点都位于地表(深度为零)，激发点间距为40 m，共51炮；接收点间距为10 m，共201个接收点；采用主频为25 Hz的雷克子波作为激发震源，频宽为10~80 Hz，炮域记录总长度为2 s，采样率为2 ms。

本次采用二维速度模型做正演测试，速度结构相对简单，输出三个SOVT就能满足成像需求。炮点正下方40°范围划为SOVT1，代表近垂直入射波场；剩余左部划为SOVT2，代表向左下入射波场；右部划为SOVT3，代表向右下入射波场(图 3)。

图 4为采用傅里叶有限差分(Fourier Finite Difference，FFD)偏移得到的三个SOVT剖面。第一水平界面照明足够，不同方向入射波场大多能被地表接收，因此三个SOVT剖面成像质量大体相当，SOVT1剖面同相轴聚焦度略好。在模型边缘位置因大角度入射波场不足，SOVT2和SOVT3剖面成像质量相对较差。第二界面水平段三个SOVT剖面成像质量与第一界面类似：靠近中心位置质量都较好，边缘位置SOVT2和SOVT3剖面成像质量较差。位于中心的凹面：SOVT1剖面对角度较小的弧面成像较好，边缘高陡弧面成像质量较差。SOVT2和SOVT3剖面分别能较好刻画高陡弧面的左半支和右半支，而对另半支弧面几乎难以正确成像。原因在于参与SOVT2剖面成像的信号主要来自弧面左半支，成像过程较少受到其他入射角度波场的影响，保证了SOVT2在凹面左半支的精确成像，SOVT3剖面也是同样道理。这三个SOVT剖面都有其相对优势的成像区域，为后续优化叠加提供了很好的数据基础。

为保证三个SOVT剖面的优化叠加，首先分析各SOVT剖面的信号聚焦特性。在第一水平界面等间距计算9个分辨率函数(范围x=200~1800 m，间隔200 m)，第二界面在水平区域同样间隔为200 m，计算6个分辨率函数(范围x=200~600 m，1400~1800 m)，在中心弧面区域增加计算密度，共计算15个分辨率函数(范围x=650~1350 m，间隔50 m)。图 5、图 6展示弧面上A点(x=700 m，z=850 m)和B点(x=1000 m，z=1000 m)分辨率函数的计算结果，并统计了相应的成像清晰度(统计范围x=0~2000 m，采样间隔1 m)和权叠加系数(表 1)。分辨率分析结果与SOVT剖面成像结果(图 4)一致：在A点，SOVT1剖面分辨率函数聚焦最好，基本无旁瓣干扰，SOVT2与SOVT3剖面聚焦结果类似，主瓣能量占比相对较小，旁瓣较大；在B点，SOVT1剖面分辨率函数聚焦最好，SOVT2剖面虽主瓣幅值较大，但旁瓣干扰也较大；总体上两者聚焦程度差异不大。SOVT3剖面分辨率函数极弱，证明对B点成像几乎无帮助，对应的权叠加系数也接近于零。

图 7左为图 4三个SOVT等权叠加剖面，等同于常规叠前FFD偏移结果，右边为优化叠加的结果。

从图 7可看出本文方法的剖面整体分辨率优于常规方法，说明优化叠加一定程度上校正了因数据频带限制等因素导致的成像误差。在弧面边缘的高陡区域，常规偏移剖面同相轴能量较弱，显然受到照明不足和偏移噪声的影响。而本文方法通过优化叠加方式相对加大了相关SOVT剖面的权重，改善了高陡弧面成像质量，初步证明本文方法的有效性。

3 实际数据测试

为进一步证明本文方法对高陡构造成像的有效性，选择位于四川盆地周边山前带的三维地震数据进行测试。资料所处位置地表起伏较大，地震激发、接收条件较差，观测系统分布不均；目的层埋深较大(约8500 m)且高陡断裂及复杂褶皱带发育。目的层照明难以保证成像质量是该区最关键的技术问题，适用于测试本文方法的性能。该三维地震数据的观测系统主要参数如表 2所示。

因主要目的层埋深较大且界面高陡，本次将数据分割为内、中、外三个圈层，其对应的入射角范围分别是0°~20°、20°~40°、40°~90°。需要注意的是，SOVT划分在偏移孔径范围内进行，外圈层最大角度90°只是说明它包含偏移孔径中的内、中圈层以外的所有地震信号。内圈层以近垂直传播波场为主，不再按方位角细分；中、外两个圈层需囊括不同走向高陡界面信号，故以方位角45°为单位分别划分为8个SOVT，从内到外一共编为17个SOVT(图 8)。

图 9为采用RTM得到的17个SOVT剖面。内圈层的SOVT1集中了近垂直传播的地震波场，传播距离短，地震波衰减较弱，剖面能量最强，能很好地刻画平缓界面，但基本不能刻画高陡界面。中、外两个圈层接收的波场传播距离较远，衰减严重，因此剖面能量相对较弱。另一方面，该地震数据测线(Inline方向)垂直地腹构造走向布设，平行构造走向的SOVT剖面(SOVT4、SOVT8、SOVT12及SOVT16)较少接收到回传地震信号，因此剖面成像能量弱于同圈层其他SOVT剖面。与正演模型测试结果类似，中线两侧的SOVT能很好地成像其对应倾向的高陡界面，为后续优化叠加提供了足够数据基础。

图 10为三维地震数据RTM剖面(左)与优化叠加剖面(右)效果对比，可见在相对平缓区域两种方法成像质量基本相当。在高陡构造区，本文方法所得剖面的同相轴的连续性及能量强度都优于常规RTM(图 10红圈内)。值得注意的是，常规RTM对红圈内高陡构造两翼成像质量存在一定差异，即左翼较差，右翼较好。提取左右两翼同深度(图 10左A、B点)的权叠加系数，从图 11可见：高陡构造左翼(A点)成像能量分布较集中(SOVT6、SOVT7)，等权叠加后有效信号占比较低，最终成像质量较差，而优化叠加后成像效果得到较大改善；高陡构造右翼(B点)因界面分布范围较大，有效成像能量分布范围较广(SOVT2、SOVT3、SOVT10、SOVT11)，等权叠加对最终成像结果影响相对较小，优化叠加后成像效果得到小幅改善(图 12)。

通过正演模型和实际数据的测试，充分证明了本文方法具备很强的高陡界面成像能力，并已在实际应用中取得良好效果。

4 讨论与结论

平缓构造区的地震反射信号在地表均匀分布，其单炮偏移道集的绝大多数炮检距、方位角段都是信号占主导，全偏数据参与叠加更有利于提高成像信噪比。在地下构造复杂，界面产状复杂多样，界面反射信号具有明显分带特征时，本文方法可通过舍弃部分偏移成果信息的方式突出有效信号，改善成像质量；而在构造相对平缓区域，常规方法与本文方法的成像效果相当。业界广泛采用的炮检距矢量片(Offset Vector Tile，OVT)数据分选方式也可分离出不同传播方向的地震波场，结合共聚焦分析也可完成优化叠加。OVT域数据多采用克希霍夫算法完成偏移，运算效率优于本文的炮域偏移算法(RTM、高斯束等)；但OVT域数据无法采用具有更高成像精度的RTM算法，可见本文方法具有更强的复杂构造刻画能力。

本文首先从炮域深度偏移出发，将每个单炮深偏道集划分为多个SOVT并叠加得到相应SOVT剖面；然后通过共聚焦理论分析各个SOVT剖面的信号聚焦特征，提取成像清晰度属性，作为后续优化叠加的权系数，最终完成各SOVT剖面的优化叠加，改善高陡界面成像质量。通过正演模型和实际地震数据测试证明本文方法能突显有效成像能量，并在一定程度上校正因观测系统不均、频带限制等导致的成像误差，最终测试结果证明本文方法具备改善原有偏移算法成像质量的能力。

总结本文方法的优势及应用局限性，有以下三点值得注意：

(1) 本文方法在常规炮域深度偏移方法上进行拓展，可同目前业界主流深度偏移算法(RTM，高斯束等)结合，只需额外增加少量运算(成像清晰度属性体可通过线性插值得到，优化叠加为叠后运算)，具备工业化推广的价值；

(2) SOVT划分方案对最后成像起关键作用，实施过程应根据面临的地震地质条件进行调整：地下构造变化剧烈区域，应适当增加SOVT数量，尽可能保证有效分离出不同倾向、倾角界面的深度偏移数据；

(3) 下一步可尝试将射线追踪算法、Poynting矢量提取等方法与SOVT划分相结合，提高对极度复杂勘探目标(高速体侵入、盐丘入侵等)地震波场的刻画精度，进一步扩展本文方法的适用范围。