0 引 言

近几十年来，地震偏移成像技术经历了射线偏移、叠后偏移、叠前偏移、逆时偏移等发展阶段，到目前为止，对偏移成像方法的理论研究已经相当成熟，但基于射线或波动方程的偏移成像理论不能用于真振幅成像、速度反演等课题，这就需要从地震反演成像的角度重新审视偏移方法(陈云峰等，2014).最小二乘偏移(LSM)是一种基于线性化反演理论的真振幅成像方法，其思路是在宏观背景模型的基础上估计出一个最优化的扰动部分对偏移结果进行迭代更新.它把成像问题当作一个反问题来处理，通过比较由偏移剖面所产生的合成数据与实际采集数据之间的相关性来判定成像结果是否准确.通过多次自动修正成像结果来提升相关性，以寻求更接近于真实的地下反射系数，从而更好地进行岩性储层成像和储层反演.

将成像问题视为最小二乘意义下的反演问题，通过最速下降法迭代使误差函数达到最小，得到横向分辨率更高的反演成像问题.其反问题可以在数据空间或者模型空间求解(Kuhl H，2003，Clapp M L，2005，Yu J H，2006)，模型空间方法不需要显式求解Hessian矩阵，但是如果不采取合理的预条件手段，迭代的收敛速度非常慢，而数据空间需要显式求解Hessian矩阵，然后对偏移结果应用Hessian矩阵的逆即可得到反演结果，相对于模型空间方法具有更快的收敛速度，但对内存的要求更高(黄建平等，2013a).进行最小二乘偏移需要求解两项：第一项为为误差泛函对模型参数的梯度，这一过程与常规偏移成像密切相关；第二项为Hessian算子，在无限孔径，无限频带和速度缓慢的假设下，Hessian矩阵是一个对角阵，但实际地震数据不可能满足这一假设，当Hessian矩阵非对角元素非零时，矩阵十分庞大，对矩阵求逆并用于反问题就变得十分困难.因此，求解一个较为精确的Hessian算子对于最小二乘偏移非常重要，自从最小二乘偏移方法提出以来，Hessian算子的计算求逆及优化一直是该方法研究的难点所在(Symes W W，2008；王华忠，2009；Xu S and Wang D，2012；任浩然等，2013).

本文通过对比分析前人的研究成果，阐述了最小二乘偏移的基本原理，指出最小二乘偏移与常规偏移的本质区别；介绍了其发展历程及研究现状；分析了最小二乘偏移实现过程中的核心问题Hessian算子及约束条件；探讨了最小二乘偏移存在的问题及今后的的发展趋势.

1 最小二乘偏移原理

基于地震波场的线性Born近似，地震检波器记录到的地震信号等效于格林函数和地震子波的卷积，可表示为

(1)

D(r，t|s，0)为观测系统记录到的地震数据，m(x)为地球反射率模型，W(t)代表震源子波，G(r，t|x，0)和G(x，t|s，0)分别代表震源(s)到散射点(x)以及散射点(x)到接收点(r)的格林函数.其中坐标系如图 1所示，(1)式矩阵表达式为

(2)

其中L为线性正演矩阵算子，它与采集形式、震源子波和速度密度模型有关，表达了地震波场传播的系统.

Kirchhoff偏移采用线性正演算子的转置来求解地下的地球反射率模型(目前常规偏移用正演算子的伴随代替其逆算子)，可表示为

(3)

由上式可知，常规偏移用正传播算子L的转置L^T来代替逆算子L^-1，当地震数据受激发和接收点采样函数以及波场的非连续性采样函数影响时，对L的近似和简化势必会降低偏移成像的分辨率及振幅的保幅性.

一方面，可以通过近似的反演公式补偿由于几何扩算引起的损失，同时对每个共深度点道集的射线条数进行归一化处理，可消除非连续采样的影响.

另一方面，在常规偏移过程中，地震记录的能量将按照椭圆轨迹外推到地下模型网格点.此时，当该点覆盖次数较高时，通过带权偏移叠加，成像位置能够较好的收敛到产生它们的地下位置点；而当覆盖次数较低时，外推能量将不能完全收敛，会产生较多的人为噪声和偏移假象(李振春，2011).

对于不规则观测系统所采集的地震数据，Kirrchoff偏移将在成像剖面上产生偏移假象，为了在最大程度上减少偏移假象，得到较为真实的地球物理模型m，引入反演的思想，构建最小二乘意义下的误差泛函为

(4)

通过对目标函数的求解，得到最小二乘解为

(5)

这种求解方法即为最小二乘偏移算法，m为最小二乘偏移结果.H=L^TL定义为Hessian矩阵(也称为照明算子、去模糊化算子等，它受波场的非连续采样函数控制).若传播矩阵L是酉矩阵(L^TL=I)或者L^T=L^－1(实际不相等)，则

(6)

(6)式所示退化为常规偏移方法.

由以上可知，常规偏移用正传播算子L的转置L^T来代替逆算子(L^TL)^－1L^T，避免了算子求逆的过程，提高了求解过程的稳定性，但是受观测系统和上覆介质的影响，缺少(L^TL)^－1项的常规偏移会导致成像分辨率低，振幅保幅性差等缺点.

本质上，常规偏移结果乘以Hessian矩阵的逆算子(L^TL)^-1，可实现一次迭代的最小二乘偏移，即

(7)

可以看出，常规偏移相当于最小二乘偏移经过Hessian矩阵滤波的结果，实际上就是模糊了的，振幅畸变了的最小二乘偏移结果.

在成像处理时，由于Hessian矩阵计算的复杂性，且直接求逆不可取，因此基于最小二乘偏移的基本原理，计算实现上主要基于最速下降法或共轭梯度迭代方法求解.其中共轭梯度法计算步骤如下：

首先，计算基于第n步迭代结果m_n的正演数据与真实数据的残差，如下式：

(8)

其次，利用偏移算子对数据残差做偏移，获得偏移成像修正量δm为

(9)

在上述给出最小二乘偏移修正量的基础上，计算出模型迭代步长表达式为

(10)

再给出第N+1步的修正的最小二乘偏移结果为

(11)

通过上述迭代过程，可以实现对成像结果的不断更新和修正，当成像误差减小至期望值或迭代次数达到给定值时，可终止迭代，获得优化的成像结果.

2 国内外研究概况

最小二乘偏移的理论基础最早由Tarantola提出(1984)，但由于计算能力的限制，导致此方法发展缓慢.最小二乘偏移反演思想由LeBras和Clayton(1988)等人提出，Lambare等(1992)进行了补充和完善.早期最小二乘偏移的方法主要基于射线理论，Nemeth等(1999)、Chavent等人(1999)在缺失数据成像、Kirchhoff真振幅成像等方面做了相关的研究.Nemeth等(1999)用惠更斯原理解释了即使数据不完整最小二乘偏移仍能得到较好的偏移效果的成因.Tamas等(1999)提出基于Kirchhoff的最小二乘偏移方法，用于偏移不规则的反射地震数据.LeBras和Clayton(1988)、Lambara等(1999)对最小二乘Kirchhoff偏移进行了不同的探讨，Nemeth(1999)对最小二乘Kirchhoff偏移方法进行了系统的研究和总结，Duquet等(2000)提出最小二乘偏移处理起伏地表照明合约和由不规则粗采样的地震波场引起的成像误差时比Kirchhoff偏移具有更大的优势.黄建平等(2012)实现了基于Kichhoff的最小二乘叠前深度偏移算法，测试了最小二乘偏移方法对复杂构造模型的偏移成像能力.黄建平等(2013)基于最小二乘Kirchhoff偏移研究了碳酸盐岩裂缝型储层地质模型的偏移方法，对不同模型在不同主频情况下进行偏移成像试算，为我国南方海相及西部探区碳酸盐岩储层的勘探开发提供了一种较为有效的针对复杂构造的成像方法.

最小二乘Kirchhoff偏移能减少偏移假象，但Kirchhoff波场传播算子的计算精度较低，不能满足实际生产的需求.Kuechl和Sacchi(1999，2001，2003)提出将最小二乘偏移应用到波场传播算子中.随后的研究主要基于单程波理论，研究重点主要集中在Hessian矩阵的快速求取、成像分辨率的提高、正则化约束、计算效率的提高等方面.Rickett(2003)通过傅里叶有限差分单程波最小二乘偏移来探讨地下照明问题.贾晓峰等(2005)推导了基于波动方程的最小二乘偏移算法.杨其强等(2008)进一步研究最小二乘偏移算法，提出了最小二乘傅里叶有限差分偏移方法，并给出了模型试算的成像结果.沈雄君等(2012)详细描述了裂步法最小二乘偏移的原理，实现算法，并通过模型进行了偏移效果测试.袁义明等(2013)提出了使用迭代再加权最小二乘的寻优方法，通过再加权的方法将稀疏优化问题转化为一系列权不断变化的线性方程组，在保证精度的前提下增加了求解的速度，对寻优过程中的线性方程组求解采用针对性的预条件化，并在权重的计算中引入伸缩变换，优化了整个寻优过程的稳定性和求解速度.黄建平等(2014)实现了一种基于分频编码方式的最小二乘裂步偏移方法.通过模型试算对比认识到：分频复用的编码方式能够使合成炮记录中的相干项两两正交，从而较好地压制超道集编码中的串扰噪声，提高成像质量；通过编码的超道集偏移在压制串扰噪声的基础上可以极大的提高成像分辨率，将分频的编码方式引入最小二乘偏移方法中，通过不断地迭代以优化成像结果，二者的结合在兼顾计算效率的同时可以得到高品质的成像剖面；采用多次、动态的编码方式迭代可以得到更好的收敛效果，引入分频策略后，成像收敛速度显著提高(相位编码方法的本质是在产生串扰噪声的相干项中使编码函数两两正交，从而达到消除相干项的目的，分频复用的编码方式对炮集进行编码，能够真正满足编码函数的正交性).在当前的计算能力下，权衡计算效率与计算精度，周华敏等(2014)提出一种基于照明补偿的单程波最小二乘偏移方法，首先利用单程波方程的稳定Borm近似广义屏波场传播算子构建反射地震数据与地下反射率间的线性算子，然后再应用线性最优化方法求解最小二乘偏移所对应的线性反问题.陈云峰等(2014)通过对波长传播算子，最小二乘成像以及Hessian算子的矩阵表达形式的详细分析，阐述了最小二乘偏移的数学物理含义，进而指出最小二乘偏移与常规偏移之间的关系，验证了最小二乘偏移较常规偏移具有更高的分辨率和保幅特性.

近期随着计算能力的大幅提高，发展了基于双程波波动方程的最小二乘逆时偏移方法(LSRTM).逆时偏移成像方法是解决复杂构造、尤其是高陡构造成像的较好方法，由于逆时偏移存在低频噪声在一定程度上影响对近地表的高精度成像，最小二乘偏移在构造恢复及能量保幅性方面具有一定优势，Dai等(2010)通过相位编码技术有效的降低了最小平方逆时偏移的计算成本并将其应用到三维数据.Yao等(2012)展示了一种最小二乘逆时偏移的矩阵实现方式，结合最小二乘公式和逆时偏移公式可对反射率进行估计，获得类似于反褶积成像条件的效果，但该方法对内存和存储量需求较高，目前不适用于三维情况.郭振波等(2014)、李振春等(2014)，黄建平等(2014)在逆时偏移成像算法的基础上，引入反演的思想实现了最小二乘逆时偏移，验证了此方法高精度保幅成像的优势.郭书娟等(2015)通过对误差泛函建立、逆时反偏移数据重构算法、Hessian逆预算条件梯度计算以及基于高斯-牛顿法的反演迭代更新方法等关键技术研究，实现了迭代最小二乘逆时偏移成像，探索了建立面向实际资料的最小二乘逆时偏移实现流程.黄建平等(2015)实现了基于静态平面波编码的最小二乘逆时偏移理论方法及处理流程，经测试对比表明，基于时移编码的平面波最小二乘逆时偏移能有效抑制低频成像噪声和串扰噪声，补偿深部成像能量，是一种有效的保幅成像策略.

3 Hessian矩阵的研究与进展

最小二乘偏移是反演框架下的偏移结果，它是对偏移成像过程和结果的校正，这种校正的桥梁就是Hessian矩阵(陈云峰等，2014).数学上，Hessian算子是误差泛函对模型参数的二阶导数，反映了误差泛函对模型变化的二次型特征；物理上，它包含了能被波动方程控制的、影响偏移振幅的所有信息.

Hessian矩阵的任意一点H_ij的求取，即为目标函数F(m)相对于该处介质参数的二阶导数，表示为

(12)

其简写为

(13)

即H由线性近似相H_a和非线性项H_n组成.Hessian逆算子在反演成像中起着重要的作用，在单参数反演中，该算子可以对照明不足位置处的参数重新聚焦，在多参数反演中，Hessian算子的非对角块矩阵中的元素描述了不同类别参数间的耦合效应(Opertoetal，2013)，使用Hessian逆算子可以减弱这种耦合效应对反演结果的影响.但由于巨大的计算量和存储量，显式计算和存储Hessian及其逆算子是不现实的，这使得解反演问题的牛顿法和高斯牛顿法得到广泛应用(Virieux and Operto，2009).在全牛顿法反演中，精确的Hessian算子包括两项：线性近似相H_a和非线性项H_n.大部分情况下，由于H_n本身计算量很大，并且在局部线性近似的条件下趋向于0，因此，首先需要研究近似的H_a(任浩然等，2013).但是，但就H_a项而言，Hessian矩阵仍然是一个巨大的矩阵，对它的求逆是最小二乘偏移中的一个难点.

基于线性Hessian矩阵的高斯-牛顿法方程最小二乘解可以表示为

(14)

假设初始反演参数m⁽⁰⁾=0，则可以得到第一次迭代过程公式为

(15)

这样，就建立了反演参数与正演算子的关系.事实上，一次迭代实现最小二乘偏移成像的过程中高斯-牛顿反演的梯度项为L^Td.将地震传播算子的共轭作用到地震数据上，即对地震数据的反传播.因此，此梯度项就是直接偏移成像的结果.可以记为m^1/2=L^Td，可以把这个介于零次迭代和一次迭代的中间结果称为是1/2次迭代结果.这个作用过程就是偏移成像，把L^T称为是偏移算子，它是正演算子L的共轭转置.m⁽¹⁾=(L^TL)^－1Ld的迭代过程就是最小二乘偏移成像.最小二乘偏移成像为真振幅成像打开了一扇大门，Hessian矩阵的逆作为加权因子作用到成像结果m^1/2上，可以看到Hessian矩阵包含了振幅修正的意义.

最小二乘偏移的难点在于求解一个精确的Hessian并求逆.在实际应用中，最小二乘函数的全Hessian矩阵非常庞大，而且难以计算，许多学者将其近似为对角矩阵，在高频近似假设记忆孔径无限大的情况下，Hessian矩阵大多为对角矩阵.在弱照明区，如盐下区域，Hessian矩阵就不再是对角阵，甚至不是对角占优矩阵.

Hessian矩阵每个对角线上的元素与介质模型的离散点相对应，其数学物理意义可以理解为：先对特定炮检对与特定绕射点形成的格林函数求自相关，然后对所有炮检对的自相关遍历求和.这个特定炮检对与特定绕射点形成的格林函数的自相关也可以理解为特定炮检对对特定绕射点的地震波“照明能量”.Hessian矩阵非对角线上的元素可以理解为：某一特定炮检对与两个不同绕射点(这两个点在模型中的位置是任意的)形成的格林函数进行互相关，然后对所有炮检对的互相关遍历求和.两个不同绕射点之间的互相关也可以理解为相邻点之间的“照明”贡献.即相邻点相距越远，其自相关函数值越大，相邻点之间“照明”贡献也越大，反之亦然(陈云峰等，2014).当一个Hessian矩阵作用于最小二乘成像结果上，至少模糊了最小二乘成像结果(实际上也改变结果的振幅值)，因此常规叠前偏移成像实际上就是模糊了的、振幅畸变了的最小二乘成像结果.由于Hessian矩阵过于庞大，远远超过目前计算机的承受能力，因此一般仅用Hessian矩阵对角元素对常规成像结果进行照明补偿，这就是常规最小二乘方法的基本实现方法.

将Hessian矩阵近似为对角阵对常规成像结果提高分辨率的能力有限.据此Albertin(2004)，Valenciano(2008)提出，计算Hessian矩阵非对角线附近的有限个元素以改善反演效果.如果要精确计算Hessian矩阵的非对角线元素，其成本也很高，因此基于各种不同的假设，提出了不同的计算Hessian矩阵非对角元素的方法.基于水平层状假设，Yu等人(2006)引入了横向不变的非对角Hessian矩阵，Guitton(2004)利用一系列非稳相滤波器近似Hessian矩阵的逆.在局部平面波假设下，Gelius(2002)，Lecomte(2008)利用射线方法在局部相位空间计算Hessian矩阵.为了避免高频近似所带来的各种问题，Xie等人(2006)采用单程波方程，通过局部平面波分解法计算相位空间Hessian矩阵.

对于叠前深度偏移，最为直接的就是利用Hessian矩阵的对角特性.从前面的分析可知，线性Hessian矩阵是一个主对角占优的带状矩阵(Ren，2011)，其特征值主要由对角元素贡献.因此，可以只利用线性Hessian矩阵的主对角元素代替整个矩阵以实现不堪精确的反演.这一假设引入的误差可以通过增加迭代次数来弥补，而一次迭代的最小二乘偏移则不能完全达到逆Hessian的效果.任浩然等(2013)对地震反演成像中的Hessian算子进行了研究，提出了两种分别适用于最小二乘偏移和全波形反演的Hessian算子简化格式.平面波Hessian算子应用于最小二乘偏移能够得到相对保真的成像结果，改善了地震偏移成像的精度；地下偏移距Hessian算子应用于全波形反演能够加快反演迭代的计算效率.陈云峰等(2014)对Hessian算子的矩阵表达形式进行了详细分析.

另外，不同类别参数间的相互耦合使多参数地震反演的非线性程度显著增加，地震波速度与各向异性参数取值数量级的巨大差异也会使反演问题的性态变差，合理使用Hessian逆算子可以减弱这两类问题对反演的影响，提高多参数反演的精度，而截断牛顿法是一种可以比较精确地估计Hessian逆算子的优化方法(Metivier et al，2013)，它在每步迭代中通过线性共轭梯度法(CG)法求解牛顿方程来获取模型增量，以此来增加准确利用Hessian逆算子的信息.Hessian矩阵的性态会影响到局部迭代优化方法的性能，当Hessian矩阵的条件数很大时(病态矩阵)，会导致局部迭代优化方法收敛缓慢甚至不收敛(Nocedal and Wright，2006).王义等(2015)采用截断牛顿法在时间域进行了VTI介质的声波双参数同时反演研究，分析了Hessian矩阵的性态.

4 最小二乘偏移的约束条件

最小二乘偏移在工业界还没有得到推广应用，究其原因主要在于两个方面：一是计算量太大，如果将LSM看作是数据空间的反演问题，其计算量约为常规偏移的2N倍(N为迭代次数)，虽然成像空间反演方式可以通过面向目标策略来提高计算速度，但是当目标区域较大时，存储成本(Hessian矩阵的存储)和计算成本仍然无法接受；第二个阻碍因素是最小二乘偏移约束条件的选取，合理的约束条件不仅可以保证反演过程稳定，改善反演效果，而且可以有效提高反演成像的收敛速度，降低迭代次数，从而也节约计算成本.但是目前地球物理学家对于LSM最优的约束条件仍然没有统一的结论(刘金玉等，2013).

自从将最小二乘反演理论引入到偏移成像以后，关于最小二乘约束条件的研究就从未停止.当已知信息不足导致反演问题存在多解性时，正则化算子所描述的约束条件非常重要.合理的正则化条件可以使反演问题的解更快地向预期方向收敛(李振春，2014).目前，国内外文献中出现的约束条件主要有：

1) 阻尼约束条件：通过加入单位矩阵保证Hessian矩阵可逆，该条件有效的前提是模型符合高斯分布，但是反射系数往往是超高斯分布的，因此该条件难以满足反射系数成像要求( Tarantola A，1984).

2) 共成像点道集的聚焦性或平滑性约束，该条件是使得反演结果在共成像点道集上沿角度坐标平滑变化或聚焦在地下零偏移距(Kuehl H，2003，Valenciano A，2006).

3) 倾角约束条件，通过拾取反射波同相轴获得局部倾角信息，利用该信息对反演结果进行约束，该方法需要介入人工干预(Prucha M，2002).

4) 预测算子，利用F-X预测滤波算子在共偏移距、角度道集对反演结果进行约束(Wang J，2008)

5) 去模糊化算子，利用两次偏移结果估计出来的滤波器作为预条件算子提高反演的收敛精度(Symes W W，2008，Aoki N，2009).

6) 变换域稀疏约束，利用成像结果在小波变换域或曲波域的稀疏分布进行约束，使得反演结果变化更加连续自然(Herrmann，2009)

7) 局部倾角约束条件，在共成像点道集引入平滑算子，在共偏移距、角度道集引入平面波构造算子进行约束，通过预条件共轭梯度法使得反偏移后的数据与输入数据之间误差达到最小(刘金玉等，2013，2015).

5 结论与讨论

最小二乘偏移计算效率低，且内存占用较大；该方法需要输入高质量的波形数据，对于数据质量较低的情况，成像收敛性较慢，有时甚至难以获得较为清晰的成像结果.最小二乘偏移的计算量约为常规偏移的2N(N为迭代次数)倍，当遇到大规模的尤其是三维地震数据时，计算量就是很大难题，用基于编码的最小二乘偏移或者引入GPU加速都可以大幅提高最小二乘偏移的计算效率，从而满足实用需求.实际地下反射系数是反射角度的函数，但目前常用的反偏移都是基于叠加成像剖面进行反偏移数据模拟，这就会导致中、远偏移距模拟数据存在比较严重的由于速度误差引起的数据时差问题，研发基于道集的反偏移算法，是完善最小二乘偏移方法所需要攻关的方向之一.

线性最小二乘偏移是对常规偏移结果的矫正，这种矫正的桥梁就是Hessian矩阵，Hessian算子中包括了控制波动方程的、影响偏移振幅的所有信息，深入研究这一算子对保幅成像和地震反演具有重要的价值.最小二乘偏移对速度模型具有较强的依赖性，迭代过程容易陷入局部极小值，合理的选取初始模型和约束条件是偏移结果准确、快度收敛的关键.

共轭梯度法对求解反演问题具有良好的收敛性，但不同的迭代方向和步长表现出不同的收敛速度，研究迭代方向的新型算法，实现全局收敛，可有效提高偏移解的稳定性.最小二乘偏移与逆时偏移，射线束偏移等算法结合是实现复杂构造保幅成像的有效方法，这是最小二乘偏移成像技术应用研究中的一个重要课题.

致谢 衷心感谢审稿专家提出的宝贵修改意见和编辑部老师的帮助.