1 引言

时间推移地震（Time-lapse seismic）是近年来新兴的一项技术，该技术根据油藏岩石物理性质变化引起的地震响应变化来监测油藏流体的变化.油藏变化主要表现在地震波速度、波阻抗、振幅和频率等方面.采用时间推移地震资料，可以直接根据不同时期的地震差异变化确定剩余油的分布^[1].探测和发现死油区，进而指导油田开发井位的合理布置，最终提高采收率^[2].

时间推移地震资料处理的关键在于如何用接收到的数据微弱变化来识别出储层的微弱变化.在时间推移地震正演的方法上，刁顺等^[3]针对于三维三相黑油模型，使用声波方程对于叠后地震记录，通过褶积模型人工合成，进行了正演模拟，并通过正演模拟进行了可行性研究.Robertsson等^{[4, 5]}提出了有限差分注入法.

在时间推移地震反演方面，Abubakar等^[6]，Gluck等^[7]分别探讨了时间推移地震波阻抗反演的可行性，Herawati^[8]计算了P波波阻抗，李来林等^[9]从最基本的地震波阻抗和地震反射系数关系式出发，推导出采用一次地震结果及测井资料计算近似时移波阻抗的新方法.李景叶^[10]等分析了利用时间推移地震AVO反演区分油藏参数变化、实现油藏定量解释的可行性，并进一步研究了多波时间推移地震AVO反演问题.陈小宏等^[11]提出了基于岩石物理模型和混合优化算法的时移地震反演压力、饱和度变化方法.此外，在时移资料互均衡方面，金龙等^{[12, 13]}也作了大量的工作.

时间推移地震（简称时移地震）不同于传统的静态地震勘探，一般而言，一次时移地震至少包含两次或两次以上的对同一区块的地震勘探.因此，一个时移地震的反演算法至少应该具备处理各个不同时期地震勘探问题的能力.第一次勘探，也称为前期勘探，是一个没有任何背景信息的地震勘探问题，因此需要大范围收敛的反演算法来求解.第二次勘探甚至包括以后的多次勘探，也称为后期勘探，是一个具有良好背景信息的地震勘探问题，需要精细、快速的反演算法来求解.鉴于此，本文首先针对传统Tikhonov正则化方法反演非连续参数时过渡光滑这一弊端，引入全变分正则化（Total Variation Regularization，简称TV正则化方法）方法.接着为了提高计算效率，引入了多尺度思想，进而构造了多尺度全变分法.而后，针对时移地震的特点，分别构造了N型多尺度全变分法和V型多尺度全变分法，以适应各个不同时期勘探的反演算法需要.最后，给出了数值模拟的结果及相关的比较和讨论.

2 全变分正则化

近几年，有一种用于图像恢复、去噪等的新方法---基于全变分的正则化方法，它可以有效地重构带“角点”的非光滑图像，即待求解可以是不连续的.这一特点特别适于用来获得波动反问题在非连续函数类中的稳定数值解.

1992年Rudin等^[14]首次引入全变分正则化方法，其主要的目的是处理一类抗干扰的图像重构问题，在处理不连续解问题及抗噪音问题时效果明显.Meyer^[15]对全变分正则化方法及其解的性质有详细的讨论.

一般全变分正则化函数可以定义成

基于对偶形式的更严格的定义是

有界变分函数空间BV（Ω）定义成BV（Ω）:={u∈L¹（Ω）：|u|_TV＜∞}.从定义可以看出空间BV（Ω）不仅包含连续函数而且包含不连续函数.

反演问题一般可以转化成求解如下泛函

(1)

的极小问题.

反问题都是不适定问题.为了数值求解的稳定性，不直接通过（1）达到极小来求解V，而是利用正则化方法在泛函（1）中引入光滑泛函，如传统的Tikhonov正则化方法就是用如下泛函

来代替目标泛函.但是一般的地震勘探反演问题，参数分布都是间断的，此时Tikhonov正则化方法^[16]对于分段（分片）常值函数的边界识别效果并不明显.为此，引入全变分正则化方法，采用

(2)

来代替（1）中的目标泛函.而（2）式中的V不要求连续，只要它是有界变差函数即可.这样，对于分段（分片）常值函数的边界识别效果有了显著的提高.

而VΔ在V=0处不可微.为了避免不可微性带来的困难，本文采用

(3)

代替（2）.其中α，β为正则系数.

易知V满足Euler方程

(4)

为此，考虑用逐步线性化的方法求解（4）.

设V^*的第k个近似值V^k已求出，为求出下一个近似V^k+1，用线性函数

代替A（V），用目标函数J_r（V）=‖L_k（V）-F‖²+代替J（V）.极小问题的相应Euler方程为：

从中解出V，令其为V^k+1，则

(5)

若非线性算子A（V）是由离散过程形成的，则（5）可以写成

(6)

其中，

3 多尺度全变分方法

多尺度算法的思想是将反问题在不同尺度上进行分解，进而求解不同尺度上的反演问题，并通过尺度之间的变化，将上一尺度求得的解作为下一尺度上反演问题求解的初值.该方法的优点在于对非连续介质模型的介质边缘具有良好的识别能力，同时计算效率也有所提高.任政勇等^[17]使用局部加密非结构化网格给出了三维电阻率计算的有限元方法，黄超等^[18]建立了可变网格与局部时间步长的交错网格的高阶差分格式，并用于弹性波模拟.田小波等^[19]构造了隐式差分多重网格算法，并将其应用于弹性波场的数值模拟中.Lu等^[20]使用代数多重网格方法计算了三维直流电阻率正演问题.

本文主要基于多重网格算法思想来构造多尺度全变分方法.多重网格是用一系列离散化和简单迭代过程来减少误差，其方法的核心是引入一系列网格G^k，在细网格上进行高频分量的光滑，而在粗网格上进行低频分量的校正.细网格上的残差通过限制算子转移到粗网格上，而粗网格上残差方程的近似解通过延拓算子插值到细网格上.借助套迭代技术由细网格到粗网格，再由粗网格依次将解返回、组合，最后在细网格上得到原问题的解.

在本文中，通过对时间推移地震所包含的不同时期勘探反问题的判断，来自适应地选择V型尺度变化或者N型尺度变化，分别构造三重网格迭代求解问题，分别为N型三重网格和V型三重网格，现分别给出两种网格类型的迭代步骤：

记三个尺度上的网格分别为Ω（h_x^l+1，h_z^l+1），Ω（h_x^l，h_z^l）和Ω（h_x^l-1，h_z^l-1），其中Ω（h_x^l-1，h_z^l-1）为粗网格，Ω（h_x^l，h_z^l）为中间网格，Ω（h_x^l+1，h_z^l+1）为最细网格.h_x^l-1=2h_x^l=4h_x^l+1，h_z^l-1=2h_z^l=4h_z^l+1.

对于N型三重网格，主要针对前期勘探反演问题，初值可以较远离真值，并且初始模型可不同于真实模型.其迭代过程如图 1所示，其中箭头代表迭代过程，虚线代表网格.

假设初值在三重网格上分别离散为V₀^l-1，V₀^l和V₀^l+1，具体N型多尺度全变分迭代算法为：

Ⅰ）k=1循环到3

（1）在粗网格Ω（h_x^l-1，h_z^l-1）上，对V_k-1^l-1使用（6）式迭代r_k1次，得到V_k^l-1；

（2）在粗网格Ω（h_x^l-1，h_z^l-1）上，计算残量，并且将残量延拓到细网格上：e^l=I_h^h_l-1l e^l-1；

（3）在细网格Ω（h_x^l，h_z^l）上，计算校正近似，并且使用（6）式迭代r_k2次，得到V_k+1^l；

（4）在细网格Ω（h_x^l，h_z^l）上，把V_k+1^l限制到粗网格Ω（h_x^l-1，h_z^l-1）上：V_k+1^l-1=I_h^h_ll-1V_k+1^l；I_h^h_l-1l

循环结束.

Ⅱ）k=1循环到2

（5）在粗网格Ω（h_x^l-1，h_z^l-1）上，对V_k+3^l-1使用（6）迭代r_k1次，得到V_k+3^l-1；

（6）在粗网格Ω（h_x^l-1，h_z^l-1）上，计算残量，并且将残量延拓到细网格上：e^l=I_h^h_ll+1I_h^h_l-1le^l-1；

（7）在细网格Ω（h_x^l+1，h_z^l+1）上，计算校正近似，并且使用（6）迭代r_k2次，得到V_k+4^l-1；

（8）在细网格Ω（h_x^l+1，h_z^l+1）上，把V_k+4^l+1限制到粗网格Ω（h_x^l-1，h_z^l-1）上：V_k+4^l-1=I_h^h_l-1lI_h^h_ll+1V_k+4^l-1；

循环结束.

Ⅲ）k=3

（9）在粗网格Ω（h_x^l-1，h_z^l-1）上，对V_k+3^l-1使用（6）迭代r_k1次，得到V_k+3^l-1；

（10）在粗网格Ω（h_x^l-1，h_z^l-1）上，计算残量，并且将残量延拓到细网格上：e^l=I_h^h_ll+1I_h^h_l-1le^l-1；

（11）在细网格Ω（h_x^l+1，h_z^l+1）上，计算校正近似，并且使用（6）迭代r_k2次，得到V_k+4^l+1；

结束.

针对后期勘探反演问题，本文构造了V型三重网格，因为可以选择前一次反演结果作为初始模型，由于油气藏开采所引起的物性变化较小，因此初始模型需较接近于真实模型，此时对算法精度要求较高.其迭代过程如图 2所示，其中箭头代表迭代过程，虚线代表网格.

假设初值在三重网格上分别离散为V₀^l-1，V₀^l和V₀^l+1，具体V型多尺度全变分迭代算法为：

Ⅰ）k=1

（1）在细网格Ω（h_x^l，h_z^l）上，计算校正近似，并且使用（6）迭代r_k2次，得到V_k+1^l；

（2）在细网格Ω（h_x^l，h_z^l）上，把V_k+1^l限制到粗网格Ω（h_x^l-1，h_z^l-1）上：V_k+1^l-1=I_h^h_l-1lV_k+1^l；

Ⅱ）k=2循环到3

（1）在粗网格Ω（h_x^l-1，h_z^l-1）上，对V_k-1^l-1使用（6）迭代r_k1次，得到V_k^l-1；

（2）在粗网格Ω（h_x^l-1，h_z^l-1）上，计算残量，并且将残量延拓到细网格上：e_l=I_h^h_l-1le^l-1；

（3）在细网格Ω（h_x^l，h_z^l）上，计算校正近似，并且使用（6）迭代r_k2次，得到V_k+1^l；

（4）在细网格Ω（h_x^l，h_z^l）上，把V_k+1^l限制到粗网格Ω（h_x^l-1，h_z^l-1）上：V_k+1^l-1=I_h^h_l-1lV_k+1^l；

循环结束.

Ⅲ）k=1循环到2

（5）在粗网格Ω（h_x^l-1，h_z^l-1）上，对V_k+3^l-1使用（6）迭代r_k1次，得到V_k+3^l-1；

（6）在粗网格Ω（h_x^l-1，h_z^l-1）上，计算残量，并且将残量延拓到细网格上：el=I_h^h_ll+1I_h^h_l-1le^l-1；

（7）在细网格Ω（h_x^l+1，h_z^l+1）上，计算校正近似，并且使用（6）迭代r_k2次，得到V_k+4^l-1；

（8）在细网格Ω（h_x^l+1，h_z^l+1）上，把V_k+4^l-1限制到粗网格Ω（h_x^l-1，h_z^l-1）上：V_k+4^l-1=I_h^h_l-1lI_h^h_ll+1V_k+4^l-1；

循环结束.

Ⅳ）k=3

（9）在粗网格Ω（h_x^l-1，h_z^l-1）上，对V_k+3^l-1使用（6）迭代r_k1次，得到V_k+3^l-1；

（10）在粗网格Ω（h_x^l-1，h_z^l-1）上，计算残量，并且将残量延拓到细网格上：el=I_h^h_ll+1I_h^h_l-1le^l-1；

（11）在细网格Ω（h_x^l+1，h_z^l+1）上，计算校正近似，并且使用（6）迭代r_k2次，得到V_k+4^l-1；

结束.

其中完全权算子I_h^h_l+1l，I_h^h_l-1l为九点限制：

双线性插值算子I_h^h_ll+1，I_h^h_l-1l为九点延拓：

4 时移地震反演模型

三维时间推移地震就是在不同的时间进行传统的二维地震勘探.描述地震波传播的二维声波方程的模型为

(7)

其中u（x，z，t）为位移函数，v（x，z）为介质在（x，z）点的速度.s（x，z，t）是震源函数，且满足s（x，z，t）=0，t＜0.x，z分别为水平方向和垂直方向的坐标.L，H为求解区域S的边界.

如果加上附加条件

(8)

这样，（7）~（8）就构成了确定速度v的声波方程反问题.

在非均匀介质中，一般求不出v的解析解，而必须用数值的方法来求解它.在这里，采用差分来代替导数进行差分离散，将反演问题转化成求解泛函的极小问题.其中是由已知的地震记录进行排列所组成向量.A是由（8）的离散方程定义的一个非线性向量值函数A：V→y^[21].

5 数值模拟

对于一个包含两次反演的时间推移地震反演问题来说，多尺度全变分法将自适应地选择所使用的网格类型.对于前期勘探反演问题，先验信息较少，所以设定的初始模型就可能偏离真实模型，此时，可以使用N型多尺度全变分方法去进行反演，得到真实模型的良好近似.对于后期勘探反演问题，可采用前期反演结果作为待反演模型的一个良好近似.因此可以使用V型多尺度全变分方法进行反演.而对于这两种算法的选择问题，由时移地震本身给出选择.

算例1：针对时间推移地震前期勘探问题，假设真实模型如图 3a所示，模拟参数为：v₁=1800 m·s^-1，v₂=2100 m·s^-1，v₃=2300 m·s^-1，采用初值模型如图 3b所示，参数为：v₁=1750 m·s^-1，v₂=2000 m·s^-1，v₃=2250 m·s^-1，图 3c是地表记录，使用该记录作为反演的附加条件.

本文分别使用N型多尺度全变分方法、基于Tikhonov正则化的正则-Gauss-Newton-多尺度方法、以及没有采用多尺度技术的单一尺度上的全变分正则化方法来进行反演.

使用的模拟参数为：L=800m，H=800m，T=1s，α=1×10^-4，β=1×10^-2，其中α，β为正则化参数，最细网格节点为32×32.分别使用三种方法进行反演，得到的结果如图 4（a，b，c）所示.为了比较三种算法的精度，定义相对误差为：

(9)

其中v为反演结果，v^*为真值.三种方法整体反演误差如图 4（d，e，f）所示，右边的标尺代表相对误差.

定义最大相对误差为：

(10)

表 1给出了三种方法的最大相对误差.从中不难看出，N型多尺度全变分方法比正则-Gauss-Newton-多尺度方法在精度上面有了很大的提高，尤其是对不连续的模型反演效果很好，是一种适用于时间推移问题初次反演的快速精细的反演算法.同时由于多重网格方法的引入，使得N型多尺度全变分方法在保持全变分正则化方法高精度的同时，大大减少了计算量.

算例2：针对时间推移地震第二次勘探问题，采用前次反演结果作为初值，这样比较接近真实值.实际上这是一个具有较好初值的反演问题，分别使用V型多尺度全变分方法，正则-Gauss-Newton-多尺度方法和单一尺度的全变分正则化方法来进行反演.

假设随着开采的进行，油藏模型从图 3a所示，变化到了图 5a所示，参数为：v₁=1800 m·s^-1，v₂=2200 m·s^-1，v₃=2000 m·s^-1，v₄=2300 m·s^-1.对应的地表记录也由图 3c所示，变化到了5c所示.

选取例一中的初次反演的真实三层斜面模型作为二次反演的初值，如图 5b所示，参数为：v₁=1800 m·s^-1，v₂=2100 m·s^-1，v₃=2300 m·s^-1.

反演使用的模拟参数为：L=800m，H=800m，T=1s，α=1×10^-4，β=1×10^-2，其中α，β为正则化参数，最细网格节点为32×32.三种方法的反演结果如图 5（a，b，c）所示.仍然采用前面（9）定义的相对误差与（10）定义的最大相对误差来比较算法精度，三种方法整体反演误差如图 5（d，e，f）所示，右边的标尺代表相对误差.

表 2给出了三种方法的最大相对误差.对比图 5（d，e，f）以及表 2，不难看出，V型多尺度全变分方法的精度是三种方法中最高的.同时，他的计算效率与正则-Gauss-Newton-多尺度方法持平，都比单一尺度的全变分正则化方法要高的多.因此，该方法是一种适用于时间推移问题二次勘探反演的快速精细的反演算法.

6 结论

本文主要构造了一种快速有效的反演方法---多尺度全变分方法.该方法首先引入全变分正则化来代替传统的Tikhonov正则化.全变分正则化方法不仅能保证反问题求解的稳定性，而且提高了对介质非连续分布区域成像的分辨能力，具有良好的边缘识别性.

接着将反问题在不同尺度上进行分解，通过对时间推移地震所包含的不同时期勘探反问题的判断，来自适应地选择V型尺度变化或者N型尺度变化.进而求解不同尺度上的反演问题，并通过尺度之间的变化，将上一尺度求得的解作为下一尺度上反演问题求解的初值.该方法的优点在于对非连续介质模型的介质边缘具有良好的识别能力，同时计算效率也有所提高.

在数值模拟中，本文分别比较了多尺度-Tikhonov正则化法，单一尺度全变分法，以及本文所构造的多尺度全变分法.在精度方面，由于全变分正则化的引入，增强了边缘识别精度，使得单一尺度全变分法与全变分多重网格法相比较，多尺度-Tikhonov正则化法具有较高的精度.在计算效率方面，多尺度上的反问题求解与单一尺度上的反问题求解相比，减少了计算量，因此多尺度-Tikhonov正则化法与多尺度全变分法具有较高的计算效率.针对两个时间推移地震反演问题的数值模拟结果也恰恰证明了这两点.

总之，本文所提出的全变分多重网格法是一种稳定、快速、精细的反演方法.该方法在时间推移地震反演中的应用进一步表明了其具有广泛的适用性.