2012, Vol. 55
2. 哈尔滨工业大学数学系,哈尔滨 150001
2. DeparLmenL of MaLhemaLics , Harbin InsLiLuLe of Technology , Harbin 150001,China
地震波形反演技术在认识地球内部结构和进行油气资源勘探研究领域具有重要的理论和实际意义.其目的是获取一个预测地震记录与实测地震记录拟合最佳的地质模型.就地震波速反演而言，线性波形反演的前提是速度场满足$V\left( {x,z} \right) = {V_0} + \delta V\left( {x,z} \right)\left( {其中\left. {\delta V\left( {x,z} \right)} \right. ＜＜{V_0}} \right)$，而这在实践中是很难被验证和满足的.事实上，地震数据记录是由地球内部物质组成和结构满足的非线性波动方程得到的，这就决定了此类问题是非线性的.1984年Albert Tarantola[1]基于声波方程和广义最小二乘理论，提出的梯度法非线性反演算法对此后近30年的地震反演理论的发展产生了深远的影响.Mora[2]延续了Albert Tarantola的工作，在迭代求梯度时提出了预条件共轭梯度算法，并利用由合成的数据反演的结果与用VSP 数据做反演得到的结果进行了对比.结果表明:在正常的地震勘探条件下，采用完全非线性地震波形反演方法，地震速度波场的所有波长分量都是可观测的，而这是线性波形反演方法无法实现的.
随着地震正演数值模拟技术和现代高性能计算机的发展与升级，在声波近似假设下，完全非线性地震波形反演方法已经被用于实际数据的处理[3-7]，实践进一步验证了非线性反演要比线性反演更接近真实.然而，由于受地震数据有限带宽和噪声的影响，加之目标函数中大量局部极值的存在，致使地震波形反演问题具有严重的不适定性.改善这种不适定性的直接方法就是增加已知信息(如假定某口井上的速度是已知的)，可以利用已知测井数据约束反演过程，提高反演过程的稳定性.文献[8-11]就是利用波动方程计算波形残差和Fréchet导数(波形残差关于速度的导数)，以测井数据为约束，通过迭代反演使合成地震记录资料与观测地震资料尽可能逼近，具有很好的精细描述储层的能力.
同伦摄动法[12-14]是求解非线性微分、积分方程的有力工具.该方法有效地结合了传统的同伦概念和摄动技术，将一个非线性问题的求解转换为无穷多个线性问题的求解，具有计算简单、收敛快速、结果稳定等诸多优势.近年来，同伦摄动法被引入到反问题的求解过程中，已有的理论和算例[15-19]表明同伦摄动法在求解各类反演问题中有巨大的应用潜力.本文针对测井约束地震波形反演问题，采用同伦摄动法迭代求解，构造了新的迭代算法.通过模型试验并与经典的Landweber迭代法进行比较，表明本文的方法反演效果更好而且算法稳定、实用.
2 测井约束地震波形反演的同伦摄动法 2.1 问题描述通常，非线性波形反演问题可以归结为如下的算子方程:
|
(1) |
其中算子F表示从参数空间X到地震观测数据空间犢的非线性映射，并假定空间X和犢为可分的Hilbert空间.实际中，我们只能得到准确数据y的含噪近似yδ，满足‖y-yδ‖ ≤δ(δ 为误差水平参数).
从测量数据应与真实数据相吻合的角度上来说，该算子方程在数学上可以转化成求解如下泛函
|
(2) |
的极小问题.
求解该问题的非线性优化方法很多，如基于导数的最速下降法[20]、共轭梯度法[21]等;基于非导数的神经网络方法[22]、遗传算法[23](GA)、模拟退火方法[24](SA)、小波多尺度方法[25]等.不同反演方法结果差别较大，不同人使用同一种方法结果也有差别，往往导致错误的地质解释.究其原因，除了地震资料品质、垂直入射假设与实际有误差等因素外，还有另外两个方面的主要影响因素:其一，由于地震数据具有带限性质，直接反演只能获得参数模型中的中频成分，缺少应有的低频和高频成分.其二，在Hardmad意义下，算子方程(1)是一个严重的不适定问题，主要表现为反演结果的不稳定性和多解性.目前，解决带限问题最有效的手段是约束反演(如测井数据的约束).而解决不适定问题主要使用正则化方法，并辅之以适当的最优化技巧[26].
2.2 测井约束模型地震资料和测井资料是地震勘探中两种最重要的资料，地震资料的分辨率较低，但它却具有范围广、横向连续性好的特点.而勘探区中的测井数据能提供井位处的地层层位的变化情况，并且具有岩性信息，但只反映地层模型坐标系中某一点的纵向变化情况.因此必须将二者结合起来，各取所长，用测井资料弥补地震资料分辨率低的缺陷[25].
速度的变化一般反映了地层岩性的变化.因此，速度参数一直是地震勘探中的一个非常重要的参数.测井约束地震波形反演就是假设某口井上的速度已知，并以此速度为约束条件，通过适当的迭代方法使地震观测数据与真实数据尽可能逼近，最终在给定误差范围内产生与地震资料特征相匹配的最优速度模型.
因此，在进行测井约束反演时，还需要给出约束点x0 处各深度的波场速度v(x0，z).为此，定义向量Vi= (vi，0，vi，1，vi，2，…，vi，n-1)T，${{\bar V}_{{i_0}}} = {\left( {{{\bar v}_{{i_0}}},0,{{\bar v}_{{i_0}}},1,{{\bar v}_{{i_0}}},2 \cdots ,{{\bar v}_{{i_0},n - 1}}} \right)^{\text{T}}}$，其中${{\bar V}_{{i_0}}}$是由位于x=x0 的声波测井资料而得到的速度. 定义容许集$\sum {{\text{ = }}\left( {V\left. {{{\bar V}_{{i_0}}} = {{\bar V}_{{i_0}}}} \right.} \right)} ,\bar V = {\left( {0,0,0, \cdots ,0,{{\bar v}_{{i_0}}},0,{{\bar v}_{{i_0},1}}, \cdots ,{{\bar v}_{{i_0},n - 1}},0,0,0, \cdots ,0} \right)^{\text{T}}}$，
|
如果$V \in \sum {} $ ，有‖DV-${\bar V}$‖ =0.
此时，极小问题(2)转化为
|
(3) |
其中λ 为约束参数，亦起到正则化的作用.
2.3 数据预处理地震声波速度数据稀疏、纵向分辨率低，原始采样测量值间隔一般在20~200m 不等，测井声波速度纵向分辨率高，一般采样间隔为0.125m，为了让两组数据分辨率匹配，必须先对数据进行预处理.
处理过程为:
(1) 对地震声波速度进行三次样条插值，使采样间隔为1m;
(2) 对测井声波速度进行多点平均，使它的分辨率和地震声波速度数据一致;
(3) 按照井深将地震声波速度和测井声波速度一一对应，形成在测井约束反演中使用的已知速度模型.
2.4 同伦摄动反演方法由于测井约束地震波形反演问题数据规模巨大，本身蕴含高度的非线性和不适定性，很多常见的迭代优化方法实现其有效求解难度不小.同伦摄动方法是一种新的、一般性地求解强非线性问题的高精度方法，已被广泛应用于各类理论与实际问题.为此，本文考虑使用同伦摄动方法求解极小化问题(3).
易知，问题(3)相应的Euler方程为:
|
(4) |
为此，引入不动点同伦曲线
|
(5) |
其中p∈[0,1]为同伦参数.当p=0时，H(V，0)=V-V0 =0.当p=1时，H(V，1)=F′(V)*(F(V)-yδ)+λD*(DV-${\bar V}$)=0.当p从0连续变到1时，同伦曲线H(V，p)从初值V0 最终求出真实解Vtrue.
现假设(5)式的解可以在p点展开
|
(6A) |
再取
|
(6B) |
将(5)式在V0 点按泰勒公式展开
|
对于p幂级数的系数相等.我们有
|
由于初值V0 已知，因此有
|
由此得到求解问题(3)的两个迭代公式，分别为一阶迭代格式:
|
(7) |
和二阶迭代格式:
|
(8) |
其中，迭代格式(7)实际上是求解测井约束优化问题(3)的Landweber迭代法.迭代格式(8)是通过对比二阶项得到的，是完全新的迭代公式，本文中将其称之为同伦摄动反演方法，它比(7)式应该更快(参见求解的类似工作[18]).本文的目的就是利用同伦摄动反演法来求解测井约束反演问题，并与Landweber迭代法(7)式进行对比.
具体计算时，若考虑地震数据噪声，迭代的停止准则为:对0 ≤k＜k* ，k*(δ)表示停止迭代步，满足
|
测井约束参数λ 在反演中起着至关重要的作用.显然，当λ￫0时，测井约束作用被降低，由于问题(2)的不适定性，必然会增大近似解的方差，导致反演过程的不稳定，最终可能促成反演失败.当λ￫∞时，测井约束作用被提升，解的光滑性和反演过程的稳定性被保证，但同时大大降低解的分辨率，最终产生一个无实际意义的解.
在反演过程中，假定在每次反演迭代中选取的测井约束参数λ 都是比较合理的，它即能保证解的方差和分辨率之间的最佳折衷，又能保证在尽可能短的时间内收敛到满足停止准则的解.因此，可构造如下的参数序列 { λk} 满足:
|
其中λ0＞0，0＜r＜1，k为迭代次数.但在实际反演中，为了体现测井约束的作用，一般选取较大的约束参数，因此，令λ0 在1.5~5之间变化.
4 数值结果为了验证本文方法的可行性，我们对两个速度模型(四层水平层状介质模型和逆冲断层带模型)，从迭代次数、相对误差、计算时间等方面进行了分析，并与经典的Landweber迭代法相对比.在正演计算中，采用位于自由表面的点脉冲震源，应用声波方程五点差分格式产生理论数据集.由于实际地震数据不可避免地带有误差，因此我们对本文中的人工合成地震数据都加入一个随机白噪声.
四层水平层状介质模型(模型Ⅰ)如图 1a所示，假设检波器分布表面在一条长2500 m 的直线上，道间距200m.本实验中，取52×26的网格，水平和垂直采样间隔均为100m，时间采样率4ms.添加噪声水平分别为0.01 和0.05，含噪合成地震记录如图 1(d，e).当测井位于水平方向第2 个网格(即i0=2)处，参数λ0=2，初速度均匀，均为2500 m/s时，同伦摄动反演方法反演得到的速度模型如图 1(b，c)所示.利用Landweber迭代法的对比实验结果如表 1所示，其中k* 表示迭代次数，τ 为CPU 运行时间，RE 为相对误差$\frac{{{V_{{k_*}}} - {V^{{\text{true}}}}}}{{{V^{{\text{true}}}}}}$
|
图 1 四层水平层状介质速度模型(模型工） (a)真实速度模型；（b)噪声水平0.01时速度反演结果；（c)噪声水平0.05时速度反演结果；(d)噪声水平0.01时地震记录；（）噪声水平0.05时地震记录. Fig. 1 The velocity of four horizontal layersmedium (model 工 (a)The true velocity model； (b) The velocity inversion result with the noisy level 0.01； (c ) The velocity inversion result with the noisy kvel 0. 05； (d) The seismograms with the noisy kvel 0. 01 ； (e) The seismograms with the noisy kvel 0. 05. |
|
表 1 模型的迭代结果 Table 1 The iteration results of models |
逆冲断层带模型(The thrust belt model，模型Ⅱ)如图 2a所示，该模型是一个复杂的典型速度模型.网格剖分、检波器分布、测井位置及λ0 取值与模型Ⅰ相同，添加噪声水平分别为0.01和0.05，含噪合成地震记录如图 2(b，c).初速度均匀取为速度模型Ⅱ的平均值，反演得到的速度模型如图 2(d，e)所示.利用Landweber迭代法对比实验结果如表 1所示.
|
图 2 逆冲断层带模型(模型n) (a)真实速度模型；（b)噪声水平0.01时速度反演结果；（c)噪声水平0.05时速度反演结果；(d)噪声水平0.01时地震记录；（）噪声水平0.05时地震记录. Fig. 2 The thrust belt model (model H ) (a) The true velocity model； (b) The velocity inversion result with the noisy level 0. 01 ； (c) The velocity inversion result with the noisy level 0. 05; (d) The seismograms with the noisy level 0. 01; (e) The seismograms with the noisy level 0. 05. |
测井约束地震反演技术已经广泛地应用于油气勘探领域，其核心是有效的迭代反演算法.随着计算机技术和数学物理等理论的不断发展，新方法和新技术不断应用于测井约束反演中，促进该技术的进一步发展.本文将地震数据与测井资料有机融合，提出了基于测井约束的同伦摄动反演方法，确保反演结果具有相对较高的纵向和横向分辨率.
通过模型试算和与经典Landweber迭代法的对比实验，可以得到如下结论:(1) 测井约束地震波形反演方法是目前提高常规地震波形反演分辨率的有效方法，波动方程反演过程是一个整体逼近的过程，据此反演得到的参数保持了测井资料的特性，具有宽频带、高分辨率的优点，并对断层的刻画具有很好的适应能力.(2) 测井约束地震波形反演方法是一种基于模型的反演技术，其关键在于模型参数的迭代修正，这也是本文的工作重点.本文将同伦摄动方法引入到反演过程中，给出了测井数据的预处理过程和测井参数的选取方法.仿真实验表明:该方法在初始模型较差的情况下，也可以获得较好的反演结果，而且具有收敛快速和计算精度高的优点.(3) 本文方法必须有测井约束资料，否则无法建立约束模型，但地下横向变化较大时，测井位置的确定会影响反演结果的精度，使方法适用于各种参数模型是未来进一步改善的方向.
| [1] | Tarantola A. Inversion of seismic reflection data in the acoustic approximation. Geophysics , 1984, 49(8): 1259-1266. DOI:10.1190/1.1441754 |
| [2] | Mora P. Elastic wave-field inversion of reflection and transmission data. Geophysics , 1988, 53(6): 750-759. DOI:10.1190/1.1442510 |
| [3] | Pratt R G, Shipp R M. Seismic waveform inversion in the frequency domain, Part II: Fault delineation in sediments using crosshole data. Geophysics , 1999, 64(3): 902-914. DOI:10.1190/1.1444598 |
| [4] | Hicks G J, Pratt R G. Reflection waveform inversion using local descent methods: Estimating attenuation and velocity over a gas-sand deposit. Geophysics , 2001, 66(2): 598-612. DOI:10.1190/1.1444951 |
| [5] | Ravaut C, Operto S, Improta L, Virieux J, et al. Multiscale imaging of complex structures from multifold wide-aperture seismic data by frequency-domain full-waveform tomography: Application to a thrust belt. Geophy. J. Int. , 2004, 159(3): 1032-1056. DOI:10.1111/gji.2004.159.issue-3 |
| [6] | Gao F, Levander A R, Pratt R G, et al. Waveform tomography at a groundwater contamination site: Vsp-surface data set. Geophysics , 2006, 71(1): H1-H11. DOI:10.1190/1.2159049 |
| [7] | Bleibinhaus F, Hole J A, Ryberg T, et al. Structure of the California Coast ranges and San Andreas fault at SAFOD from seismic waveform inversion and reflection imaging. Journal of Geophysical Research , 2007, 112: B06315. DOI:10.1029/2006JB004611 |
| [8] | 冯国峰, 韩波, 刘家琦. 二维波动方程约束反演的大范围收敛广义脉冲谱方法. 地球物理学报 , 2003, 46(2): 265–270. Feng G F, Han B, Liu J Q. Widely convergent generalized pulse-spectrum methods for 2-D wave Equation Inversion. Chinese J. Geophys (in Chinese) , 2003, 46(2): 265-270. |
| [9] | 傅红笋, 韩波. 二维波动方程速度的正则化-同伦-测井约束反演. 地球物理学报 , 2005, 48(6): 1441–1448. Fu H S, Han B. A regularization homotopy method for the inverse problem of 2-D wave equation and well log constraint inversion. Chinese J. Geophys. (in Chinese) , 2005, 48(6): 1441-1448. |
| [10] | Han B, Fu H S, Liu L. A homotopy method for well log constraint waveform Inversion. Geophysics , 2007, 72(1): R1-R7. |
| [11] | 王守东. 井间地震资料全变差正则化波形反演. 石油物探 , 2011, 50(3): 234–240. Wang S D. Total variation regularization waveform inversion of crosswell seismic data. Geophysical Prospecting for Petroleum (in Chinese) , 2011, 50(3): 234-240. |
| [12] | He J H. Homotopy perturbation technique. Comput. Meth. Appl. Mech. Eng. , 1999, 178(3-4): 257-262. DOI:10.1016/S0045-7825(99)00018-3 |
| [13] | He J H. A coupling method of a homotopy technique and a perturbation technique for non-linear problems. Int. J. Non-Linear Mech , 2000, 35(1): 37-43. DOI:10.1016/S0020-7462(98)00085-7 |
| [14] | Sharma P R, Methi G. Applications of homotopy perturbation method to partial differential equations. Asian Journal of Mathematic and Statistics , 2011, 4(3): 140-150. DOI:10.3923/ajms.2011.140.150 |
| [15] | Shakeri F, Dehghan M. Inverse problem of diffusion equation by He's homotopy perturbation method. Phys. Scr. , 2007, 75(4): 551-556. DOI:10.1088/0031-8949/75/4/031 |
| [16] | Ganji D D, Sadighi A. Application of He's homotopy-perturbation method to nonlinear coupled systems of reaction-diffusion equations. Int. J. Nonlinear Sci. Numer. Simul. , 2006, 7(4): 411-418. |
| [17] | Javidi M, Golbabai A. A numerical solution for solving system of Fredholm integral equations by using homotopy perturbation method. Appl. Math. Comput. , 2007, 189(2): 1921-1928. |
| [18] | Abdulaziz O, Hashim I, Chowdhury M S H. Solving variational problems by homotopy-perturbation method. Int. J. Numer. Meth. Eng. , 2008, 75(6): 709-721. DOI:10.1002/nme.v75:6 |
| [19] | Cao L, Han B, Wang W. Homotopy perturbation method for nonlinear ill-posed operator equations. Int. J. Non. Sci. Numer. Simul. , 2009, 10(10): 1319-1322. |
| [20] | 崔岩, 王彦飞, 杨长春. 带先验知识的波阻抗反演正则化方法研究. 地球物理学报 , 2009, 52(8): 2135–2141. Cui Y, Wang Y F, Yang C C. Regularizing method with a priori knowledge for seismic impedance inversion. Chinese J. Geophys. (in Chinese) , 2009, 52(8): 2135-2141. |
| [21] | 许琨, 王妙月. 声波方程频率域有限元参数反演. 地球物理学报 , 2001, 44(6): 852–864. Kun X, Wang M Y. Finite element inversion of the coefficients of acoustic equation in frequency domain. Chinese J. Geophys. (in Chinese) , 2001, 44(6): 852-864. |
| [22] | 周辉, 何樵登, 徐世浙. 人工神经网络非线性地震波形反演. 石油物探 , 1997, 36(1): 61–70. Zhou H, He Q D, Xu S Z. Nonlinear seismic waveform inversion using the artifiual neural network. GPP (in Chinese) , 1997, 36(1): 61-70. |
| [23] | Sambridge M, Drijkoningen G. Genetic algorithms in seismic waveform inversion. Geophy. J. Int. , 1992, 109(2): 323-342. DOI:10.1111/gji.1992.109.issue-2 |
| [24] | 姚姚. 地球物理非线性反演模拟退火法的改进. 地球物理学报 , 1995, 38(5): 643–650. Yao Y. Improvement on nonlinear geophysical inversion simulated annealing. Chinese J. Geophys. (in Chinese) , 1995, 38(5): 643-650. |
| [25] | 马坚伟, 朱亚平, 杨慧珠. 二维地震波形小波多尺度反演. 工程数学学报 , 2004, 21(3): 109–113. Ma J W, Zhu Y P, Yang H Z. 2-dimension seismic waveform wavelet multiscale inversion. J. Engi. Math. (in Chinese) , 2004, 21(3): 109-113. |
| [26] | 王彦飞. 反演问题的计算方法及其应用. 北京: 高等教育出版, 2007 . Wang Y F. Computational Methods for Inverse Problems and Their Applications (in Chinese). Beijing: Higher Education Press, 2007 . |