应用反射理论的最小二乘逆时偏移成像

段心标 王华忠 邓光校

(①同济大学海洋与地球科学学院波现象与智能反演成像研究组，上海 200092；②中国石化石油物探技术研究院，江苏南京 211103； ③中国石化西北油田分公司，新疆乌鲁木齐 830011)

0 引言

从地震数据中准确提取并利用复杂的地下岩性与储层信息，是油气勘探的关键环节。地震波成像是连接地震观测数据与储层解释的桥梁，常规偏移成像技术[1-5]主要目的是获得能精细描述地质构造的地震成果资料，用于部署和落实井位。

随着油气勘探程度的日趋深入和地震成像技术的不断进步，直接估算地下弹性参数的反演成像方法受到越来越多的关注。早在20世纪60年代，Backus等[6-7]从数学角度研究了地球物理反问题，奠定了早期地球物理反演的理论基础。Tarantola[8-9]基于Bayes估计理论，建立了一套地震波全波形反演的理论方法体系，推动了地震反演成像技术的快速发展。然而，全波形反演具有很强的非线性特征，难以由观测数据直接反演得到较准确的介质弹性参数，无法满足油藏描述的精度需求。实际中应用的依然是经典的地震波成像处理和储层描述流程，而如何获得更高偏移成像质量是其中重要环节。

最小二乘偏移[10-18]是常规偏移的后续发展，其核心思想是在线性反演理论下求解模型空间的精确解，相当于在常规偏移构造成像基础上，考虑了振幅补偿和波形校正，从而得到更保幅和更高分辨率的成像结果。

经典最小二乘偏移以散射波线性化表达为基础，认为模型参数是速度(或慢度)的高波数扰动，其成像实质是估计介质的散射强度。Zhang等[19-20]提出基于逆时偏移叠加剖面的反偏移方法，并形成相应的最小二乘偏移技术，但并未给出反偏移方法的详细理论依据。陈生昌等[21-22]提出地震反射波数据偏移成像方法和相应的反射波动方程最小二乘偏移方法，并引入“反射率”概念，即最小二乘偏移的实质是估计反射面上的反射率。本文认为：因地下介质以层状为主，获取高保真反射系数是地震成像走向储层刻画的关键，估算不同角度的反射系数是最小二乘偏移的目标。然而，考虑到反演不同角度反射系数的计算量巨大，且高质量偏移叠加成像依然是当今实际应用的主要目标，因此以最小二乘偏移获得忽略角度信息的平均反射系数具有现实意义。

本文首先简介线性最小二乘偏移方法基础理论，重点讨论散射理论线性化表达方法和反射理论线性化表达方法，并分别建立散射强度、反射系数与地震数据的线性关系式； 然后从偏移成像实际应用需求出发，选择基于反射理论并忽略角度信息的Kirchhoff正演反偏移表达方法，实现与实际应用需求相一致的最小二乘逆时偏移技术； 最后，通过模型数据和实际资料验证本文方法的效果。

1 方法原理

1.1 线性最小二乘偏移方法

基于Bayes估计理论，地震反演成像可表示为后验概率密度最大化问题； 在地震反演通常假设条件下，进而可转化为如下目标函数最小化问题

(1)

式中：m为地下介质参数，如速度、散射强度、反射系数或波阻抗等；dobs为观测数据；L(·)表示波场传播算子。该式非线性反演问题求解难度大，需对非线性反演问题逐步进行线性化，依次反演参数的低、中、高波数信息。最小二乘偏移是一种线性化的地震反演成像方法，用于估计介质参数的高波数成分，其目标函数可表示为

(2)

式中L为矩阵形式的线性化波场传播算子。

令目标函数梯度为零，则有如下的法方程

LTLm=LTdobs

(3)

最小二乘偏移成像结果表示为

(4)

式中：LTL为Hessian矩阵，可用H表示；LTdobs是对观测数据进行常规偏移成像。最小二乘偏移是Hessian矩阵的逆H-1作用于常规偏移成像，能够消除Hessian矩阵导致的成像模糊化效应，提高成像分辨率和保幅性，减少成像假象。

由于“目标函数梯度为零”仅是一种理想假设，且Hessian矩阵非常庞大，其计算和存储都很困难，Hessian矩阵的求逆更难以完成，因而无法直接求解法方程以实现最小二乘偏移。为此，通常采用迭代法求解最小二乘偏移问题，假设目标函数在初始值附近满足二次型，从初始解逐步逼近真实解。

1.2 地震波场的线性化表达

1.2.1 散射理论的Born正演线性化表达

一阶Born近似假设下，散射波场可写成

(5)

(6)

进而，散射波场可表示为

G0(x′,xs,ω)dx′

(7)

式(7)定义了一个Fredholm第一类积分方程，它描述了一个线性问题。其积分核可写成

(8)

在散射假设下，可认为地表接收到的携带地下信息的信号就是散射波场记录，则可将式(7)写成如下线性算子表达式

Lm=dobs

(9)

此式中的m设定为慢度平方扰动，则该式表达了散射强度与地震数据之间的线性关系。

1.2.2 反射理论的Kirchoff正演线性化表达

反射波理论假设下，地下界面Ω上入射波场引起的二次场可表示为

(10)

地下界面Ω上波场的法向导数为

(11)

利用WKBJ近似方法[23]，有

(12)

由Green定理，反射面上方地面上任意一点x的反射波场可写成如下积分式

(13)

式中Gr(x,x′,ω)是反射波传播的Green函数(图1)。

图1 波场传播示意图

将式(10)～式(12)代入式(13)，可得

(14)

根据程函方程关系，式(14)可进一步写成

Gs(x′,xs,ω)Gr(x,x′,ω)dΩ

(15)

式中s(x′)为慢度矢量。

式(15)建立了地震数据与不同方向反射系数之间的联系[24-25]。然而，在缓变速度结构中，不同角度的反射波能量较均衡； 在速度剧变结构中，不同角度的反射系数难以估计。同时，考虑到直接基于式(15)进行波场模拟的计算量巨大，因此建立地震数据与平均反射系数之间的关系就具有现实意义。若不考虑背景速度对振幅的影响，式(15)就退化为

(16)

式中：D表示散射体；R′表示平均反射系数。

在反射假设下，认为地表接收到的携带地下信息的信号就是反射波场记录。式(16)同样可写成线性算子(式(9))表达形式，此时m表示忽略地下角度信息的平均反射系数，其积分核为

(17)

1.3 最小二乘逆时偏移实现方法

逆时偏移技术现已在实际处理中得到广泛应用，其成像基础是Claerbout[26]给出的成像方法：波场外推+成像条件。成像条件为：成像点处，入射波的到达时等于反射波的出发时。显然，常规逆时偏移是估计地下介质的反射系数。由于地下介质以层状为主，以定位反射界面的位置或产生角度反射系数道集为目标的反射系数成像是当前地震成像的主流需求。

如前所述，由于角度反射系数反演的计算量巨大，因而实际应用中，最小二乘偏移的目的是在常规逆时偏移基础上进一步提高偏移叠加成像质量，使成像结果分辨率更高，成像假象更少。

在反演成像方法中，反偏移采用Kirchhoff正演思路，具体迭代反演步骤如下。

(1)以设定的偏移速度场对观测数据dobs作逆时偏移成像，作为最小二乘偏移的初始值m0。

(2)对第i次迭代后偏移成像值mi(i=0,1,2,…)做反偏移，正演模拟预测数据dcal，反偏移公式为

(18)

(3)计算观测数据dobs与预测数据dcal的残差波场δd，当残差满足设定条件时迭代停止。

(4)对δd进行逆时偏移成像，得到目标泛函的梯度gi。

(5)计算共轭梯度方向

(19)

(6)计算迭代步长

(20)

(7)更新偏移成像值

mk+1=mk-μk+1hk

(21)

(8)转到步骤(2)，展开新一轮迭代。

由于迭代算法计算量巨大，三维最小二乘逆时偏移中采用了GPU加速计算技术，在此不再赘述。

2 数值模型试验及实际应用

首先利用二维Sigsbee2a模型(图2a)做数值试验，验证本文方法效果； 进而将三维最小二乘逆时偏移技术应用于实际地震资料成像处理。

图2 Sigsbee2a模型及其成像结果

2.1 Sigsbee2a模型数据

为了便于计算，把Sigsbee2a经典模型速度值单位转为m/s，并重新定义速度网格，新网格的尺寸为3200×1200，间隔为5m×5m。对该模型正演模拟120炮数据，炮间距为100m，检波点间隔为20m，炮检距范围是-2000～2000m，记录长度为7s。

利用本文所提最小二乘逆时偏移方法进行反演迭代，得到第20次迭代结果(图2c)，并与常规逆时偏移剖面(图2b)进行对比； 图3是目标泛函迭代误差分析曲线，可见反演逐步收敛。从整体剖面及盐丘左下局部放大图(图4)的对比可见，相对于常规逆时偏移方法，本文方法的成像效果(图4b)显著改善，分辨率变高，深部成像更清晰。

图3 迭代误差曲线

图4 Sigsbee2a模型的常规逆时偏移(a)和最小二乘逆时偏移(b)成像的局部放大

2.2 实际资料

本次试验所用实际数据取自中国西北M油田，测试处理的目的是提高储层的成像精度，力求对奥陶系缝洞体及小断裂精确成像。所选数据满覆盖面积为401km2，共48000炮，成像面元是30m×30m。利用100GPU节点(每节点4块K10卡)3次迭代，共用时59天。

图5和图6是两条不同测线的常规逆时偏移与最小二乘逆时偏移成像剖面对比，可见最小二乘逆时偏移剖面中“串珠”成像收敛更好，对深层大断裂、层间小断层的刻画更清晰，T78(奥陶系蓬莱坝组顶界面)以下成像能量更强，成像分辨率得到提高。对比奥陶系目的层的成像频谱(图7)，可见最小二乘逆时偏移成像拓宽了频谱(图7b，频带展宽5Hz)； T74(一间房组顶界面)下40～120ms均方根振幅，也证实本文方法成像结果(图8b)对“串珠”刻画更清晰。

以上结果均展示了最小二乘逆时偏移成像在提高成像分辨率、振幅均衡性及成像精度等方面的效果。

图5 L1测线的常规逆时偏移(a)和最小二乘逆时偏移(b)成像剖面

图6 L2测线的常规逆时偏移(a)和最小二乘逆时偏移(b)成像剖面

图7 常规逆时偏移(蓝线)和最小二乘逆

3 结论与讨论

最小二乘偏移是一种线性化的地震反演成像方法，线性近似可以是“背景速度+速度扰动”模式，也可以是“背景速度+反射系数”模式。由于地下介质以层状为主，现今业界常规偏移成像的主要目标是反射系数。本文分析了两种线性化模式的差异，选择基于忽略角度信息的反射系数的反偏移正问题表达方法，提出一种估计平均反射系数的最小二乘逆时偏移方法。理论模型测试结果表明本方法能提高成像分辨率，深部成像更清晰。将研发的三维最小二乘逆时偏移技术应用于中国西部探区资料处理，所得“串珠”成像比常规逆时偏移结果收敛更好，深层大断裂、层间小断层的刻画更清晰，成像分辨率获得提高，充分体现了最小二乘逆时偏移成像在提高成像分辨率和振幅均衡性及提高成像精度方面的优势。

图8 常规逆时偏移(a)与最小二乘逆时偏移(b)的均方根振幅平面属性对比

然而，目前最小二乘逆时偏移实际应用效果仍不甚理想，主要原因包括反偏移得到的无噪声一次波模拟数据与实际观测记录差异较大、实际地震子波未知、现有建模精度达不到最小二乘逆时偏移要求及三维技术计算量大等。为了降低最小二乘逆时偏移对速度模型和地震子波的依赖，近年也发展了动态拉伸时差校正方法和空变子波估计方法，一定程度上减弱了反演的病态性。但考虑到反问题的复杂性，这些应用瓶颈问题并未得到有效解决，它们会持续伴随最小二乘逆时偏移的发展过程，有待逐步克服和解决。