基于预条件共轭梯度法的宽带约束波阻抗反演

尹德强

中国电力工程顾问集团东北电力设计院有限公司，长春 130021

0 引言

宽带约束反演是一种地震岩性反演方法，能够充分利用原始地震资料提供地层的空间分布和断层等有效信息，利用测井资料补偿地震资料缺失的低频和高频信息，拓宽了地震资料的频带，有效解决地震数据的带限问题，恢复出宽频带波阻抗的空间分布[1]。地震数据和波阻抗参数之间呈现非线性数值关系，实质要解决非线性反演问题[2--3]。而非线性反演理论体系并不完善，反演过程中的收敛速度慢、解的稳定性差及其反演精度低等存在的问题常限制其应用[4]。

针对求解雅克比系数矩阵病态问题，研究者对共轭梯度算法改进策略不同。采用正则化加入阻尼因子降低矩阵条件数；对系数矩阵采用预处理方法降低矩阵条件数；对搜索方向的改进措施，采用P--R--P方法和DY方法等，但多数用于稀疏约束反褶积间接反演波阻抗。波阻抗反演存在带限问题，主要采取测井资料进行约束，对于横向阻抗急剧变化反演效果较差。笔者以地震资料解释的标志层作为控制点，依据测井资料提供的井点信息，以井点处进行外推内插，建立合适的宽频带初始波阻抗模型。以地震数据与合成地震记录的残差为基础建立反演目标函数，利用参数置换法对原始波阻抗进行参数变换，反演解加入硬约束条件，对转化的拟线性矩阵方程采用改进的预条件共轭梯度算法进行数值求解以获得高分辨率波阻抗模型。

1 构建反演目标函数

采用阻尼最小二乘方差拟合函数建立地震反演目标函数：地震数据与合成地震记录的残差为基础，采用模型参数最小长度解为约束：

J=s-Wr2+μrTr

(1)

式中：s为地震数据；r为反射系数序列；W为子波矩阵；μ为阻尼因子。

(GTG+μI)·ΔZ=GTΔS

(2)

针对式(2)求解波阻抗矩阵方程，可采用线性化算法进行求解。但在实际迭代过程中，系数矩阵G每一次迭代会因初始模型的改变而发生变化，而系数矩阵G往往矩阵的条件数很大，呈现病态特征。同时搜索方向不一定是极值方向，会影响收敛速度和反演解的精度。依据参数置换法，设置波阻抗对数为新的参数，L(i)=ln[(Zi)]，反射系数与波阻抗映射关系为：

r=DL

(3)

式中：r=(r1,r2,…,rn)T,L=(L1,L2,…,Ln)T，为波阻抗对数矩阵；D为常数矩阵，形成的矩阵如式(4)所示：

(4)

原始问题目标函数式(1)可转化为：

J=‖s-WDL‖2+μDLTDL

(5)

目标函数式(5)在新参数波阻抗对数进行泰勒级数展开，可得矩阵方程：

(GTG+μI)·L=GTS

(6)

式中：G为子波矩阵W与常数矩阵D形成的系数矩阵；S为地震数据矩阵；阻尼因子μ与地震数据和模型的方差矩阵有关。

由式(6)可知，在应用共轭梯度法进行迭代求解时，系数矩阵G为常数矩阵，合理的选择阻尼因子μ，能够降低系数矩阵的条件数，一定程度上减少在反演迭代过程中对收敛速度和搜索方向的影响。应用传统Fletcher--Reeves方法迭代求解式(6)矩阵方程时，算法后期迭代过程中，容易产生连续小步长，收敛速度很慢。

2 改进的预条件共轭梯度法

共轭梯度法是利用初始点处的梯度方向构造一组共轭方向，沿着共轭方向进行搜索目标函数的极值，具有收敛速度快和二次终止性。传统Fletcher--Reeves方法算法依赖于初始模型的选择，易陷入局部极值。针对系数矩阵的病态性，常用的方法有正则化加入阻尼因子以降低矩阵条件数；对系数矩阵进行不完全Cholesky分解、不完全的LU分解等降低系数矩阵条件数；采用重开始策略，令第n次迭代结果为新的初始点重开始，以使其最终达到收敛[6--7]。

采用Cauchy稀疏分布构造反射系数约束项，建立误差的最小二乘拟合函数[8--9]，矩阵方程为：

(7)

对式(7)进行极小化并将矩阵方程进行分解，预条件矩阵R以乘积的形式作用于系数矩阵为：

(8)

笔者采用一种改进的近似预条件共轭梯度算法，令：G=WTW，d=WTs，改进的措施是对负梯度d-Gr进行预条件处理，具体的迭代步骤如下所示：

采用改进的预条件共轭梯度法可改善系数矩阵的病态特征，提高收敛速度和反演精度，一定程度上可提高程序运行效率。

3 宽带约束波阻抗反演

宽带约束波阻抗反演以原始地震数据资料为基础，在先验知识的约束下，以地质信息、测井资料为约束条件，建立宽频带的初始波阻抗模型；采用合适的数值方法求解，最终获取具有宽频带的高分辨率波阻抗模型[10--12]。

采用改进的预条件共轭梯度算法，使得搜索方向靠近负梯度方向，改进的搜索方向为：

p(k+1)=g(k+1)+βkp(k)

(9)

式中：p(k+1)为新的搜索方向；g(k+1)为负梯度方向；参数βk更改为：

(10)

从系数矩阵G本身性质出发，其中∂Ji/∂Lj表示数据si在解分量Lj方向上的变化程度，则合成地震记录J在解分量Lj总的变化率可以表示为：

(11)

式中：T为衰减因子，T=a·(k-1)；a为常数；k为迭代次数。

本文所采用的宽带约束波阻抗反演主要步骤：

①由于地震资料属于带限信号，宽带约束反演的关键是建立宽频带的初始波阻抗模型，补偿缺失的低频信息；以地震资料解释的标志层作为控制点，依据测井资料提供的井点信息，以井点处进行外推内插，建立合适的宽频带初始波阻抗模型。可在地震道中选取井资料进行约束。

②在反演过程中，要对反演解进行约束，不然会造成反演解偏离初始猜测的模型处太远，获得的相对波阻抗期望特征相同；但绝对波阻抗值会出现严重误差，不是根据初始猜测波阻抗模型推出数据趋势得到。反演过程中对反演解加入硬约束：

Lmx(i)≤L(i)≤Lmz(i)

(12)

式中：Lmx(i)为约束反演解最小值；Lmz(i)为约束反演解最大值，取值的范围可约束在偏离初始猜测值±15%。

③建立反演目标函数，采用参数置换法对原始波阻抗进行参数变换，得到式(6)近似线性系统；采用更改的预处理共轭梯度法求解，对初始波阻抗模型进行迭代修改，直到求解的模型波阻抗合成的地震记录与原始地震数据最佳匹配为止，进而获得最优的高分辨率波阻抗模型，宽带约束波阻抗反演流程见图1。

图1 宽带约束波阻抗反演流程图Fig.1 Inversion flow chart of broadband constrained wave impedance

4 模型试算及结果分析

设计一个12层单道理论速度模型及一个6层初始速度模型，相关系数仅为0.785(图2)，采样时间为2 ms，密度常数1.0，采样点数235个。选取地震子波为零相位阻尼余弦子波，采样率2 ms，子波的长度为40 ms，主频为60 Hz，地震记录是由地震子波与反射系数相褶积形成，加入5%的随机噪声。

利用fortran语言编制计算机程序进行波阻抗反演，其反演效果如图3所示：

图3 地震记录中加入5%随机噪声反演结果对比Fig.3 Inversion results comparison with 5% random noise in seismic record

由图3可知，在原始地震记录中加入5%的随机噪声，反演结果与理论模型相关系数可达0.997，剖面与绝对波阻抗值均与理论值大体相同；对比采用稀疏约束反褶积法[13--14]，由迭代求出的反射系数存在误差，小的反射系数误差经过合并后，递推公式产生的波阻抗会产生很大累积误差，造成深部波阻抗值偏离真实值很大。而对于矩阵求逆法、不完全Cholesky分解及不完全的LU分解，存在迭代过度，收敛不稳定等问题。笔者直接对波阻抗模型进行反演，改进搜索方向收敛速度更快，反演精度高，抗燥性强。

以一个61道砂泥岩互层理论速度模型模拟实际油气储层。该模型大致可分为5层：地表为速度较低的第四系覆盖层；第二层为泥岩层；第三层为楔状砂岩体储层，局部含有油气资源；而第四层泥岩层局部存在逆断层，使得下部第五层含水砂岩向上侵入，各层的岩性及速度值见表1。选取的地震子波为阻尼余弦子波，地震道的长度为0.4 s，采样点数为200个，采样率为2 ms，理论速度模型如图4所示。依据第35道速度值建立虚拟井数据，构建与理论模型相关性较差的初始速度模型(图5)。

表1 理论模型速度参数Table 1 Speed parameters of theoretical model

图4 砂泥岩互层岩性油气藏模型Fig.4 Lithologic reservoir model in sand-shale interbed

图5 构建速度初始模型Fig.5 Initial model of established velocity

设各个岩层的密度均为常数，采用褶积模型合成地震记录：阻尼余弦子波与理论速度模型相褶积合成人工地震记录剖面，并在原始地震记录中加入能量比为5%的随机噪声，则形成的拟实际地震剖面如图6所示。在无噪音及5%随机噪声下，采用改进的预条件共轭梯度法进行宽带约束波阻抗反演，反演重构地下地质模型，恢复含油储层和断层等地质构造信息，并反演出较为真实的速度模型。

图6 加入5%随机噪声的合成地震记录Fig.6 Synthetic seismic record with 5% random noise

如图7所示，无噪声情况下，反演得到的速度模型较好反映出楔状含油储层和下部地层存在的逆断层位置，准确描述出各个反射层位，相关系数可达0.998，各道反演速度点的绝对误差最大仅为真实值的±4%。

图7 无噪音反演得到的速度模型Fig.7 Velocity model obtained from inversion without noise

如图8所示，存在5%的随机噪声干扰时，反演得到速度模型也能够描述出储层和逆断层的位置，由于噪音强度分布不同，造成各道反演得到的绝对速度值存在偏差，但最大偏差值仅为真实值的±5%。基于宽带约束波阻抗反演研究方法对比中，笔者依据系数矩阵的特征巧妙构造预处理矩阵，既不改变原始初始问题，又能利用预处理技巧；在迭代过程中，改进搜索方向，对反演解进行加入硬约束条件，充分利用测井信息，使得反演解是由初始猜测的波阻抗模型推出的数据趋势而得到。

图8 5%随机噪音反演得到的速度模型Fig.8 Velocity model obtained from inversion with 5% random noise

由图9可知，基于预条件的宽带约束波阻抗反演较大的偏差均集中在15～25道之间，是由于给出初始模型相关性较差的原因，以单道进行反演，没有考虑相邻道相互约束，造成相邻道数据存在不同误差。通过拟实际数据验证，采用的改进的预条件共轭梯度算法应用于宽带约束反演方法，验证反演算法的正确性和优越性。

图9 反演结果的平均误差曲线Fig.9 Average absolute error curves of inversion results

5 结论

(1)宽带约束反演建立优化的初始波阻抗模型，可以减少多解性问题，依赖于地震资料的品质和高信噪比，适用于井资料丰富的地质区域。

(2)用改进的预处理共轭梯度法，对反演解加入硬约束条件，实现基于模型反演的宽带约束波阻抗反演。在5%随机噪声干扰时，反演最大偏差值仅为5%，且真实反映储层和逆断层位置。笔者采用的反演算法收敛速度快、反演精度高、数值稳定性好，具有一定的抗噪性。

(3)线性化的数值算法进行迭代求解，反演的结果常与初始模型选择有关，通过复杂理论模型试算，偏差值较大集中在15～25道，收敛到局部极值，未考虑相邻道约束。因此，发展和完善完全非线性反演理论是地震反演方法发展必然趋势。